{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Benchmark code for 3d correlation function\n",
    "\n",
    "This code requires cluster_toolkit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pyccl as ccl\n",
    "import cluster_toolkit\n",
    "#cluster toolkit package is available at http://cluster-toolkit.readthedocs.io/en/latest/source/installation.html\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "bench_ind = 3\n",
    "\n",
    "if bench_ind == 1:\n",
    "    cosmo = ccl.Cosmology(\n",
    "        Omega_c=0.25, Omega_b=0.05, h=0.7, sigma8=0.8, \n",
    "        n_s=0.96, w0=-1.0, wa=0.0, Omega_g=0,\n",
    "        transfer_function='bbks')\n",
    "elif bench_ind == 2:\n",
    "    cosmo = ccl.Cosmology(\n",
    "        Omega_c=0.25, Omega_b=0.05, h=0.7, sigma8=0.8, \n",
    "        n_s=0.96, w0=-0.9, wa=0.0, Omega_g=0,\n",
    "        transfer_function='bbks')        \n",
    "elif bench_ind == 3:\n",
    "    cosmo = ccl.Cosmology(\n",
    "        Omega_c=0.25, Omega_b=0.05, h=0.7, sigma8=0.8, \n",
    "        n_s=0.96, w0=-0.9, wa=0.1, Omega_g=0,\n",
    "        transfer_function='bbks')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = np.logspace(np.log10(5.e-5), 3., 730000) # Wavenumber\n",
    "h = 0.7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CCL power spectrum\n",
    "pk = []\n",
    "for n in range(6):\n",
    "    pk.append(ccl.nonlin_matter_power(cosmo, k, 1./(n+1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate CCL xi and benchmark xi for r = 0.1 - 100 with 40 bins\n",
    "nr1 = 40\n",
    "r1 = np.logspace(-1, 2, nr1)\n",
    "\n",
    "xi1 = []\n",
    "for n in range(6):\n",
    "    xi1.append(ccl.correlation_3d(cosmo, 1./(n+1), r1))\n",
    "\n",
    "xi_toolkit1 = []\n",
    "for n in range(6):\n",
    "    xi_toolkit1.append(cluster_toolkit.xi.xi_mm_at_r(h*r1, k/h, pk[n]*h*h*h, exact=True))    \n",
    "        \n",
    "# calculate CCL xi and benchmark xi for r = 50 - 250 with 100 bins to check agreement in teh BAO peak region\n",
    "nr2 = 100\n",
    "r2 = np.logspace(np.log10(50), np.log10(250), nr2)\n",
    "\n",
    "xi2 = []\n",
    "for n in range(6):\n",
    "    xi2.append(ccl.correlation_3d(cosmo, 1./(n+1), r2))\n",
    "\n",
    "xi_toolkit2 = []\n",
    "for n in range(6):\n",
    "    xi_toolkit2.append(cluster_toolkit.xi.xi_mm_at_r(h*r2, k/h, pk[n]*h*h*h, exact=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# write benchmark xi to file\n",
    "with open('model%d_xi.txt' % bench_ind,'w') as f:\n",
    "\n",
    "    f.write('#  [0] r (Mpc; comoving), [1] xi(r,z=0.0), [2] xi(r,z=1.0), [3] xi(r,z=2.0), [4] xi(r,z=3.0), [5] xi(r,z=4.0), [6] xi(r,z=5.0)' + '\\n')\n",
    "\n",
    "    for i in range(140):\n",
    "        col = []\n",
    "        s = ''\n",
    "        if i < 40:\n",
    "            col.append(\"{:.18e}\".format(r1[i]).ljust(25))\n",
    "            for n in range(6):\n",
    "                col.append(\"{:.18e}\".format(xi_toolkit1[n][i]).ljust(25))\n",
    "        else:\n",
    "            col.append(\"{:.18e}\".format(r2[i-40]).ljust(25))\n",
    "            for n in range(6):\n",
    "                col.append(\"{:.18e}\".format(xi_toolkit2[n][i-40]).ljust(25))\n",
    "        s = ' '.join(col)\n",
    "        f.write(s + '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "r                        xi(CCL)                   xi(benchmark)            Delta(r^2 xi)       Delta xi / xi_benchmark\n",
      "1.00000e-01              9.13885e+02              9.16199e+02              -2.31471e-02             -2.52643e-03             \n",
      "1.19378e-01              7.68389e+02              7.69716e+02              -1.89134e-02             -1.72422e-03             \n",
      "1.42510e-01              6.44119e+02              6.43962e+02              3.18136e-03              2.43254e-04              \n",
      "1.70125e-01              5.34180e+02              5.34738e+02              -1.61577e-02             -1.04400e-03             \n",
      "2.03092e-01              4.41434e+02              4.41346e+02              3.65583e-03              2.00827e-04              \n",
      "2.42446e-01              3.61355e+02              3.60953e+02              2.35976e-02              1.11221e-03              \n",
      "2.89427e-01              2.92533e+02              2.92692e+02              -1.33039e-02             -5.42614e-04             \n",
      "3.45511e-01              2.34771e+02              2.34896e+02              -1.50246e-02             -5.35800e-04             \n",
      "4.12463e-01              1.86435e+02              1.86349e+02              1.47169e-02              4.64216e-04              \n",
      "4.92388e-01              1.45962e+02              1.45925e+02              9.05724e-03              2.56007e-04              \n",
      "5.87802e-01              1.12633e+02              1.12560e+02              2.51080e-02              6.45606e-04              \n",
      "7.01704e-01              8.54888e+01              8.54752e+01              6.72653e-03              1.59825e-04              \n",
      "8.37678e-01              6.37750e+01              6.37642e+01              7.56558e-03              1.69087e-04              \n",
      "1.00000e+00              4.66669e+01              4.66823e+01              -1.53626e-02             -3.29088e-04             \n",
      "1.19378e+00              3.34944e+01              3.35053e+01              -1.55918e-02             -3.26540e-04             \n",
      "1.42510e+00              2.35885e+01              2.35906e+01              -4.28709e-03             -8.94810e-05             \n",
      "1.70125e+00              1.63185e+01              1.63184e+01              2.01773e-04              4.27217e-06              \n",
      "2.03092e+00              1.11338e+01              1.11329e+01              3.73905e-03              8.14267e-05              \n",
      "2.42446e+00              7.54132e+00              7.54487e+00              -2.08882e-02             -4.70996e-04             \n",
      "2.89427e+00              5.12649e+00              5.12499e+00              1.25711e-02              2.92821e-04              \n",
      "3.45511e+00              3.52741e+00              3.52908e+00              -1.99483e-02             -4.73502e-04             \n",
      "4.12463e+00              2.48550e+00              2.48401e+00              2.54811e-02              6.02970e-04              \n",
      "4.92388e+00              1.79456e+00              1.79423e+00              7.81466e-03              1.79645e-04              \n",
      "5.87802e+00              1.32690e+00              1.32649e+00              1.43232e-02              3.12518e-04              \n",
      "7.01704e+00              9.97268e-01              9.97372e-01              -5.15032e-03             -1.04874e-04             \n",
      "8.37678e+00              7.57602e-01              7.57521e-01              5.65924e-03              1.06466e-04              \n",
      "1.00000e+01              5.76894e-01              5.76763e-01              1.31717e-02              2.28372e-04              \n",
      "1.19378e+01              4.36847e-01              4.36969e-01              -1.73464e-02             -2.78556e-04             \n",
      "1.42510e+01              3.27605e-01              3.27644e-01              -7.72376e-03             -1.16074e-04             \n",
      "1.70125e+01              2.41917e-01              2.41923e-01              -1.59306e-03             -2.27518e-05             \n",
      "2.03092e+01              1.75213e-01              1.75183e-01              1.22660e-02              1.69757e-04              \n",
      "2.42446e+01              1.23922e-01              1.23929e-01              -4.27461e-03             -5.86803e-05             \n",
      "2.89427e+01              8.53563e-02              8.53289e-02              2.29466e-02              3.21030e-04              \n",
      "3.45511e+01              5.69251e-02              5.69318e-02              -8.02919e-03             -1.18139e-04             \n",
      "4.12463e+01              3.66221e-02              3.66250e-02              -4.93160e-03             -7.91483e-05             \n",
      "4.92388e+01              2.25658e-02              2.25750e-02              -2.21267e-02             -4.04273e-04             \n",
      "5.87802e+01              1.32290e-02              1.32243e-02              1.63932e-02              3.58782e-04              \n",
      "7.01704e+01              7.28325e-03              7.28097e-03              1.12484e-02              3.13759e-04              \n",
      "8.37678e+01              3.70677e-03              3.70348e-03              2.30823e-02              8.88210e-04              \n",
      "1.00000e+02              1.68299e-03              1.68555e-03              -2.55455e-02             -1.51556e-03             \n",
      "5.00000e+01              2.15880e-02              2.15979e-02              -2.48889e-02             -4.60949e-04             \n",
      "5.08195e+01              2.06085e-02              2.05990e-02              2.45351e-02              4.61192e-04              \n",
      "5.16524e+01              1.96330e-02              1.96396e-02              -1.75770e-02             -3.35452e-04             \n",
      "5.24990e+01              1.87178e-02              1.87197e-02              -5.16084e-03             -1.00028e-04             \n",
      "5.33594e+01              1.78327e-02              1.78309e-02              5.22497e-03              1.02918e-04              \n",
      "5.42340e+01              1.69712e-02              1.69797e-02              -2.50872e-02             -5.02319e-04             \n",
      "5.51229e+01              1.61626e-02              1.61590e-02              1.11543e-02              2.27177e-04              \n",
      "5.60263e+01              1.53661e-02              1.53693e-02              -1.02233e-02             -2.11911e-04             \n",
      "5.69446e+01              1.46158e-02              1.46146e-02              4.06669e-03              8.58124e-05              \n",
      "5.78779e+01              1.38817e-02              1.38876e-02              -1.95608e-02             -4.20470e-04             \n",
      "5.88265e+01              1.31856e-02              1.31882e-02              -9.03699e-03             -1.98013e-04             \n",
      "5.97907e+01              1.25232e-02              1.25213e-02              6.84458e-03              1.52909e-04              \n",
      "6.07706e+01              1.18812e-02              1.18818e-02              -2.23642e-03             -5.09665e-05             \n",
      "6.17666e+01              1.12640e-02              1.12666e-02              -9.97464e-03             -2.32057e-04             \n",
      "6.27790e+01              1.06726e-02              1.06792e-02              -2.62325e-02             -6.23263e-04             \n",
      "6.38079e+01              1.01164e-02              1.01132e-02              1.28125e-02              3.11167e-04              \n",
      "6.48537e+01              9.57219e-03              9.57422e-03              -8.55573e-03             -2.12463e-04             \n",
      "6.59167e+01              9.06216e-03              9.06043e-03              7.48027e-03              1.90011e-04              \n",
      "6.69970e+01              8.56761e-03              8.56611e-03              6.74262e-03              1.75361e-04              \n",
      "6.80951e+01              8.09103e-03              8.09702e-03              -2.77823e-02             -7.39966e-04             \n",
      "6.92112e+01              7.64777e-03              7.64514e-03              1.25833e-02              3.43604e-04              \n",
      "7.03455e+01              7.21519e-03              7.21665e-03              -7.22947e-03             -2.02441e-04             \n",
      "7.14985e+01              6.80662e-03              6.80522e-03              7.12285e-03              2.04747e-04              \n",
      "7.26703e+01              6.41456e-03              6.41552e-03              -5.03768e-03             -1.48691e-04             \n",
      "7.38614e+01              6.04059e-03              6.04129e-03              -3.78703e-03             -1.14904e-04             \n",
      "7.50719e+01              5.69086e-03              5.68768e-03              1.78848e-02              5.57947e-04              \n",
      "7.63023e+01              5.35018e-03              5.34771e-03              1.43327e-02              4.60345e-04              \n",
      "7.75529e+01              5.02906e-03              5.02632e-03              1.64470e-02              5.44053e-04              \n",
      "7.88240e+01              4.72220e-03              4.72088e-03              8.22836e-03              2.80526e-04              \n",
      "8.01159e+01              4.43004e-03              4.42970e-03              2.16777e-03              7.62433e-05              \n",
      "8.14290e+01              4.15318e-03              4.15331e-03              -8.95675e-04             -3.25235e-05             \n",
      "8.27636e+01              3.89092e-03              3.88978e-03              7.78257e-03              2.92092e-04              \n",
      "8.41201e+01              3.64323e-03              3.64033e-03              2.05747e-02              7.98720e-04              \n",
      "8.54988e+01              3.40072e-03              3.40342e-03              -1.97344e-02             -7.93213e-04             \n",
      "8.69001e+01              3.17681e-03              3.18020e-03              -2.55957e-02             -1.06579e-03             \n",
      "8.83244e+01              2.96736e-03              2.96766e-03              -2.39090e-03             -1.03273e-04             \n",
      "8.97720e+01              2.76462e-03              2.76690e-03              -1.83463e-02             -8.22760e-04             \n",
      "9.12433e+01              2.57767e-03              2.57551e-03              1.79882e-02              8.38923e-04              \n",
      "9.27388e+01              2.39806e-03              2.39680e-03              1.08689e-02              5.27270e-04              \n",
      "9.42588e+01              2.22672e-03              2.22570e-03              9.09621e-03              4.59992e-04              \n",
      "9.58037e+01              2.06431e-03              2.06618e-03              -1.71361e-02             -9.03605e-04             \n",
      "9.73739e+01              1.91207e-03              1.91391e-03              -1.74184e-02             -9.59851e-04             \n",
      "9.89698e+01              1.76949e-03              1.77187e-03              -2.32701e-02             -1.34079e-03             \n",
      "1.00592e+02              1.63726e-03              1.63692e-03              3.47588e-03              2.09851e-04              \n",
      "1.02241e+02              1.51235e-03              1.51028e-03              2.16666e-02              1.37242e-03              \n",
      "1.03916e+02              1.39305e-03              1.39197e-03              1.16731e-02              7.76583e-04              \n",
      "1.05619e+02              1.28175e-03              1.28006e-03              1.88196e-02              1.31793e-03              \n",
      "1.07351e+02              1.17286e-03              1.17405e-03              -1.36501e-02             -1.00889e-03             \n",
      "1.09110e+02              1.07424e-03              1.07552e-03              -1.52666e-02             -1.19232e-03             \n",
      "1.10898e+02              9.85145e-04              9.82959e-04              2.68842e-02              2.22388e-03              \n",
      "1.12716e+02              8.95877e-04              8.96339e-04              -5.86076e-03             -5.14649e-04             \n",
      "1.14563e+02              8.14749e-04              8.15793e-04              -1.36935e-02             -1.27892e-03             \n",
      "1.16441e+02              7.40124e-04              7.39497e-04              8.50385e-03              8.48140e-04              \n",
      "1.18349e+02              6.70468e-04              6.69523e-04              1.32322e-02              1.41102e-03              \n",
      "1.20289e+02              6.03175e-04              6.02709e-04              6.74821e-03              7.73800e-04              \n",
      "1.22261e+02              5.40904e-04              5.41836e-04              -1.39403e-02             -1.72119e-03             \n",
      "1.24264e+02              4.83696e-04              4.84466e-04              -1.18858e-02             -1.58881e-03             \n",
      "1.26301e+02              4.30739e-04              4.30437e-04              4.81953e-03              7.01908e-04              \n",
      "1.28371e+02              3.81689e-04              3.81559e-04              2.14925e-03              3.41814e-04              \n",
      "1.30475e+02              3.35926e-04              3.35058e-04              1.47616e-02              2.58797e-03              \n",
      "1.32614e+02              2.92335e-04              2.92422e-04              -1.54330e-03             -3.00098e-04             \n",
      "1.34787e+02              2.52844e-04              2.53306e-04              -8.38964e-03             -1.82306e-03             \n",
      "1.36996e+02              2.17124e-04              2.16685e-04              8.23640e-03              2.02531e-03              \n",
      "1.39242e+02              1.83477e-04              1.83149e-04              6.35881e-03              1.79075e-03              \n",
      "1.41524e+02              1.52523e-04              1.52194e-04              6.59108e-03              2.16223e-03              \n",
      "1.43843e+02              1.24058e-04              1.24053e-04              1.07661e-04              4.19440e-05              \n",
      "1.46201e+02              9.79308e-05              9.74066e-05              1.12049e-02              5.38168e-03              \n",
      "1.48597e+02              7.39943e-05              7.35130e-05              1.06268e-02              6.54661e-03              \n",
      "1.51033e+02              5.21110e-05              5.17918e-05              7.28253e-03              6.16425e-03              \n",
      "1.53508e+02              3.21498e-05              3.20064e-05              3.38020e-03              4.48171e-03              \n",
      "1.56024e+02              1.39863e-05              1.39324e-05              1.31219e-03              3.86890e-03              \n",
      "1.58581e+02              -2.49755e-06             -2.94895e-06             1.13516e-02              -1.53069e-01             \n",
      "1.61180e+02              -1.74134e-05             -1.75428e-05             3.36058e-03              -7.37380e-03             \n",
      "1.63822e+02              -3.08672e-05             -3.09893e-05             3.27850e-03              -3.94202e-03             \n",
      "1.66507e+02              -4.29589e-05             -4.33776e-05             1.16084e-02              -9.65252e-03             \n",
      "1.69236e+02              -5.37832e-05             -5.40283e-05             7.02139e-03              -4.53750e-03             \n",
      "1.72010e+02              -6.34294e-05             -6.35953e-05             4.91057e-03              -2.60976e-03             \n",
      "1.74829e+02              -7.19817e-05             -7.22102e-05             6.98126e-03              -3.16307e-03             \n",
      "1.77694e+02              -7.95198e-05             -7.96801e-05             5.06145e-03              -2.01177e-03             \n",
      "1.80607e+02              -8.61185e-05             -8.62383e-05             3.90635e-03              -1.38868e-03             \n",
      "1.83567e+02              -9.18483e-05             -9.20047e-05             5.27198e-03              -1.70049e-03             \n",
      "1.86575e+02              -9.67754e-05             -9.69787e-05             7.07456e-03              -2.09563e-03             \n",
      "1.89633e+02              -1.00962e-04             -1.01407e-04             1.60031e-02              -4.38839e-03             \n",
      "1.92741e+02              -1.04467e-04             -1.04879e-04             1.52956e-02              -3.92580e-03             \n",
      "1.95900e+02              -1.07345e-04             -1.07751e-04             1.55846e-02              -3.76878e-03             \n",
      "1.99111e+02              -1.09648e-04             -1.09906e-04             1.02441e-02              -2.35104e-03             \n",
      "2.02375e+02              -1.11423e-04             -1.11864e-04             1.80562e-02              -3.94116e-03             \n",
      "2.05692e+02              -1.12716e-04             -1.13128e-04             1.74407e-02              -3.64386e-03             \n",
      "2.09063e+02              -1.13569e-04             -1.13747e-04             7.80731e-03              -1.57039e-03             \n",
      "2.12489e+02              -1.14021e-04             -1.14217e-04             8.88302e-03              -1.72248e-03             \n",
      "2.15972e+02              -1.14109e-04             -1.14313e-04             9.55732e-03              -1.79244e-03             \n",
      "2.19512e+02              -1.13867e-04             -1.14222e-04             1.71255e-02              -3.11155e-03             \n",
      "2.23109e+02              -1.13327e-04             -1.13723e-04             1.97160e-02              -3.48284e-03             \n",
      "2.26766e+02              -1.12519e-04             -1.12764e-04             1.25821e-02              -2.16984e-03             \n",
      "2.30483e+02              -1.11471e-04             -1.11873e-04             2.13729e-02              -3.59633e-03             \n",
      "2.34260e+02              -1.10208e-04             -1.10455e-04             1.35653e-02              -2.23792e-03             \n",
      "2.38100e+02              -1.08754e-04             -1.09001e-04             1.40255e-02              -2.26971e-03             \n",
      "2.42002e+02              -1.07130e-04             -1.07359e-04             1.33709e-02              -2.12660e-03             \n",
      "2.45969e+02              -1.05358e-04             -1.05753e-04             2.38664e-02              -3.73024e-03             \n",
      "2.50000e+02              -1.03457e-04             -1.03673e-04             1.34967e-02              -2.08297e-03             \n"
     ]
    }
   ],
   "source": [
    "# print some values\n",
    "n = 0 # redshift    \n",
    "print(\"r                        xi(CCL)                   xi(benchmark)            Delta(r^2 xi)       Delta xi / xi_benchmark\")\n",
    "for i in range(140):\n",
    "    col = []\n",
    "    s = ''\n",
    "    if i < 40:\n",
    "        col.append(\"{:.5e}\".format(r1[i]).ljust(25))\n",
    "        col.append(\"{:.5e}\".format(xi1[n][i]).ljust(25))\n",
    "        col.append(\"{:.5e}\".format(xi_toolkit1[n][i]).ljust(25))\n",
    "        err = r1[i]*r1[i]*(xi1[n][i]-xi_toolkit1[n][i])\n",
    "        col.append(\"{:.5e}\".format(err).ljust(25))\n",
    "        rel_diff = (xi1[n][i]-xi_toolkit1[n][i])/xi_toolkit1[n][i]\n",
    "        col.append(\"{:.5e}\".format(rel_diff).ljust(25))\n",
    "    #s = col[0] + col[1] + col[2] + col[3] + col[4]\n",
    "    else:\n",
    "        col.append(\"{:.5e}\".format(r2[i-40]).ljust(25))\n",
    "        col.append(\"{:.5e}\".format(xi2[n][i-40]).ljust(25))\n",
    "        col.append(\"{:.5e}\".format(xi_toolkit2[n][i-40]).ljust(25))\n",
    "        err = r2[i-40]*r2[i-40]*(xi2[n][i-40]-xi_toolkit2[n][i-40])\n",
    "        col.append(\"{:.5e}\".format(err).ljust(25))\n",
    "        rel_diff = (xi2[n][i-40]-xi_toolkit2[n][i-40])/xi_toolkit2[n][i-40]\n",
    "        col.append(\"{:.5e}\".format(rel_diff).ljust(25)) \n",
    "    s = col[0] + col[1] + col[2] + col[3] + col[4]\n",
    "        \n",
    "    print(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEoCAYAAACQD2yQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2debxVZdXHvwuZQUQBL4MCoRKSAgoOYCWKKZYImqGGiHOOae9bDmVlaWXZ+2ZvKWqDcw5pDqgpDlzUQAZNUMRw4qoMDigiMxfW+8c6m7Pv4cx3n3l9P5/zOefs/ey9n7vvPvu31/CsR1QVx3Ecx4mCFqXugOM4jlM9uKg4juM4keGi4jiO40SGi4rjOI4TGS4qjuM4TmS4qDiO4ziR4aLilC0icriIPCMiKiIrRGSKiLwgIv8Wke+LSMss99NKRM4WkQdF5Lpm9GeoiKwUkaH57iPH4/UWkatE5BUR6RNb1jN2LsYmtO0Wa/tnEXlXRA6JLT9SRH4kIi+KyNMi0r4YfXdqFxcVp2xR1anANbGv16vqGFU9ELgK+A1wa5b72QTcAgwEOuTSBxFpF/q6AngG+CyXfeSLqr4LLAD2AiS2eDUwDVie0PxaYI6qngHcBXQQkR2BG4FfAmOAlUBWQuw4+eKi4pQ762Lvm4MFqno/8C/g2yLSI5udqOp6YGkuBxaRScABoX0sVtVjVfXNXPbTTJaFv6jqKlU9TlVnhfrZFjiOmNip6iWq+ggwDtiixlJV/aaqripi350axEXFqVQ+jL3vUoidi8gwIKmrTES2K8QxcyGhD12B1kma9SpSdxxnK24KOxVH7IY6GNgIvBVa3gb4MVAH7As0AGep6scp9tMduBL4D7AbsAW4KPY+AXOVnSMiX8PccKfEXhcB9SLyVeABYDtgpKq+HLOc7sXcVufEDnUx8AVgb8zy+o6qbu13Qp9aA78G2mOWVc/Qunaxfp0C/Am4VUQOx6wUYn09CrgZGAt8HdhRRH4LrFPVH4uIJOsP8A5wOHAqcF9s26OArwJvJzuvmDvtdOAkYDIwHDgZeAMYp6rvx/rdEfg+5sLbC1gPnKuqn6XqT6rz41QAquovf5XtCxgJKHBF7PvOwJ8xd9hFCW2vBwbHPrcHPgLuDa2vB24Jfb8fuDv2eTvsJnlS7Hvf2HFHxr53wOISW5fFll+CiVvn0LI7gXah9UfHPgvwCjA7zd97K/C/oe8Xx47ZF2gD7BP7fkqoTd/EfsWWXwEsTliWtD9AW0xAFBPKEcBNQO9U5zW2fe/YNo8Ae8b+P+8C18bat8BiQPvFvnfCROWX+Zwff5X/yy0Vp1I4RkRGYTGOe4F9VHV+sFJEegInAqtE5MTY4tmkt8YfAIIYgwKfA/2SNVTVNSLySpJVf8Ju3hOA62LWz3JVXRezOi4DbhKREbH2i4BuItJSVRvDOxKRgdiT/sDQ4tmhPmwA/m0P97mTrj9Ao6o+G9v306o6A5iR7ryqKcG7sW3uU9WFseM8jwkMmMXUUlXnxP6GVSIyDngr1/PjVAYuKk6l8ADm1nkJc8G8k7B+EObiuTTbHarqHSKyk4hcgLm8WpNjnFFVPxGRvwNnYDGYCcSz0voBOwCXq+rGLHZ3aOz93Vz6kAPZ9md96HPO5xWz3NrEPn8VWBJeqaqPA4jIgCz741QQHqh3KgZVbQAmAgOwFOEwHYCdk43DiD0Rb4OIHAk8ij1lX0c80yxXbgCGiMi+wJ4hCypIX/5CkmO3SVwWat85z35kItf+BNvkdF4TaIW5yKLqj1PmuKg4FYWqPoaNuzhWRC4LrXoDi4ucHm4vIqdi2VEkLG8D/A2LuQRpu+HfQzDRUEZfU8xVNJ/YWJHQqrcwC+ishGOPBr6UZFdBqvJXMx0zT3LtD+R4XpPwOnCgiOyTsP0xefbHKXPc/eWUO51i7+En158ABwJXiciLqjpVVefHfPm/EZFWmN//YKCtqgbjU1phN0iwp+TOwH6xJ+4jgC5ADxHpCnyKCcsAEVmAucZaxbZN9oR+I5YhdlSwQFVXisjfgO+JyFrgcSzQPlhVz0yyj0cx19cvRGQuJjJHxNYdLCIriLumwn0IPic+3bcOL8vUHxEJRLVlaJu05zVU1SAsvmE34m3AT4EnRORnwGLgm1gyQq7nx6kESp0p4C9/pXoBo7BBjordYM8IresGvI8F2i8Ftgd2BR7D3FjvYTezFsSftNdisZhvxPbxv8AaYCE2UPD3WGbTpNj6P8f2/7+Yi+b3sb78AxiS0NcdgGuS/A2dMYtoNfAB8H9A+zR/8xexUfurgReA82LfTwP6AD+K9eFZ4MuYK/AG4llbh8b28w0srVmxFOh+6foTO39Xxto/F+wntk2683pZbJvHgD2Ar2Cj/T8FjohtPwwTo3XAXOAr+Z4ff5X/S2L/WMdxHMdpNh5TcRzHcSLDRcVxHMeJDBcVx3EcJzJcVBzHcZzIcFFxHMdxIqPmx6l07dpV+/btm9e2a9asoUOHnOZ8cvLAz3Ph8XNcHKrpPL/44osfq2q3xOU1Lyp9+/Zl7ty5eW1bX1/PyJEjo+2Qsw1+nguPn+PiUE3nWUQaki1395fjOI4TGS4qjuM4TmS4qDiO4ziRUbOiIiJjROSmzz77rNRdcRzHqRpqVlRUdYqqnrXDDjuUuiuO4zhVQ82KiuM4jhM9LiqO4zhOZLioOI7jVAqNjXDMMZDn2LpiUJWiIiJ7l7oPjuM4kfPGG/Dgg3DSSaXuSUqqTlREZH9sxjzHcRynyFSdqKjqbGxKWMdxnOpCpNQ9yEjZiIqItBURz+91HMfJRBlPA19yURGRFiIyCVgE7JewrpeITBaRc0TkNhHZqzS9dBzHKQMCS6WMRaUcqhR3BaYDu4YXiogADwOXqepUEZkOPCoiewCHAhcn2dexqrqq0B12HMcpCRXg/iq5qKjqhwCy7ck6DBgI1MfavSYim4BxqnofMLWI3XQcxykfythSKbn7Kw0HAe+o6sbQskWYlZISEdkX6CYihxeyc47jOEWnjMUkoOSWShq6A4nVHlcCu6TbSFVfAtJOrSYiZwFnAdTV1VFfX59XB1evXp33tk72+HkuPH6Oi0Nzz3P7hgb2B9auW8fsMv1/lbOoNAKbEpZFYlmp6k3ATQDDhg3TfGdiq6ZZ3MoZP8+Fx89xcWj2eV6wAID2bduW7f+rnN1fy4DEFOPOwNIodu6l7x3HqTi2bLH3MnaDlbOoTAP6iUi4jwOIBe6bi5e+dxyn4ti8udQ9yEhZiEqCcATMBBqAkbE2A7BYyZSIjumWiuM4lYVbKpkRkZ2BS2NfJ4jIngCqqsBYYJKInAdcBhylqmuiOK5bKo7jVByBqJQxJQ/Ux8ap/DL2Slz3FjAp9vW6YvbLcRyn7HD3V/ni7i/HcSoOd3+VL+7+chyn4nBRcRzHcSLD3V/li7u/HMepONxSKV/c/eU4TsVRAdlfNSsqjuM4FYe7v8oXd385jlNxuPurfHH3l+M4FYeLiuM4jhMZ7v5yHMdxIsMtlfLFYyqO41Qcnv1VvnhMxXGcisPdX47jOE5kuPvLcRzHiQwXFcdxHCcy3P3lOI7jRIZbKuWLZ385jlNxePZX+eLZX47jVBzu/nIcx3Eiw91fjuM4TmS4qDiO4ziR4e4vx3EcJzLcUnEcx3Eiw7O/yhdPKXYcp+Jw91f54inFjuNUHO7+chzHcSLDRcVxHMeJDHd/OY7jOJHhlorjOI4TGZ795TiO40SGu7+Kh4h0F5F/iEiDiPys1P1xHMeJHHd/ZYeItBWR5ub2HgKMB/YGviMinZvfM8dxnDIiEJUydoOVVFREpIWITAIWAfslrOslIpNF5BwRuU1E9sqwu/tVtVFVVwGvAesK1G3HcZzSELi/ytgN1rLEx+8KTAd2DS8UEQEeBi5T1akiMh14VET2UNXGZDtS1Y2xbbsBT6nqhsJ23XEcp8gEFoqLSnJU9UMA05AmHAYMBOpj7V4TkU3AOBFZBVycZHfHAp8DY4BfF6jLjuM4paMC3F+ltlRScRDwTmB9xFgEHKqq5wJTk20kIscC96jqZhHZVVXfK0JfHcdxioO7v/KmO5BY6XElsEuqDUTkHOBS4HMRaQ2cByQVFRE5CzgLoK6ujvr6+rw6uXr16ry3dbLHz3Ph8XNcHJp7nndraLBYQWMjS8aNY/Epp7Cpc3nlJJWrqDQCmxKWpU0qUNXJwORsdq6qNwE3AQwbNkxHjhyZRxehvr6efLd1ssfPc+Hxc1wcmn2eH3po68deDz1Er86d4ZZbmt2vKCmLlOIkLAMSU4w7A0ujOoCXvnccp+JIdHttKL98pHIVlWlAPxEJ928AscB9FHjpe8dxKo7EAH0ZDoIsuagkCEfATKABGBlrMwDoAEwpXs8cx3HKjERLpQyzwEoaUxGRnYEzYl8niMgSVV2oqioiY4GfiMiewP7AUaq6JsJjjwHG7L777lHt0nEcp7AkioiLSlNi41R+GXslrnsLmBT7el0Bjj0FmDJs2LAzo9634zhOQXBRKV/cUnEcp2JYsgSWLasI91fJYyqlwgP1juNUDMcdB/vtB4sXN13ugXrHcRwnZ15/3d6nTWu6vAxH1tesqPg4FcdxKoLGRli7Nvm6VMtLSM2Kiru/HMepCN58EzZuTL5uTWQJsZFRs6LiOI5TEXz4ob136WLv228fX7d6dfH7k4GaFRV3fzmOUxEEGV67xqadCntX3FIpH9z95ThORZAoKmFcVBzHcZycCESlTx97f//9+LrPP49/Xr8ejjgC5swpXt+S4KLiOI5TzgRjUfbbb9t1GzaYmAC89RZMnQqjRhWvb0lwUXEcxylnAkslVfWPIC4ciMvnn8Mdd8DZZ5dkcGTNiooH6h3HqQgCUdluOzjsMHvdeSdMnGjLg3tYOL4ycSLceCPMmhVfpmr7UoXzz4cZMwrS3ZoVFQ/UO45TEQSiIgJPPmmvb38bxo+35StX2nuy9OIHH7T3f/0LvvIVE6Z16+C66+DLXy5Id2u2oKTjOE5FEIhKiwQbIJibPtFSqauDDz6wz59/btuHBSQY9yJSkO7WrKXiOI5TEQRxkURRCbwsp55q1kdgqfTuHW+zdi289FLT7d55pzD9jJGTqIhIVxH5vojMFJEPRWSdiLwpIveKyNcL1UnHcZyaJZOlsmQJ3HdfclFZs2bb2MnbbxemnzGydn+JyCRgAjAb+D9gFbAe6AR0B06ItTlXVVcUoK+R4vOpOI5TEWQSFYA2beLur/AgybVr4ZNPmm4XWCoFcn9lJSoicg7QoKqHp2k2WUR6ABeLyBWqui6SHhYIn/nRcZyKIByoD9OxY/zzypVmqbRoAV27xpevXRuPuQQUWFQyur9EpAvwnKo+Flr2FxEZmthWVZcBPwL2jrSXZciTT8KMGV1K3Q3HcaqNn/4Ujj8+/j1VTCUsCitWmKXSsWPTdmvXxrPDAgrs/sooKqq6QlVfTVjcHViWon2jqs6OonPliipcdRVcccWXePLJUvfGcZyq4uc/h3vvjX9P5f4CeP55e1+xwiyVDh0yi0o5BepD/AA4MubuqjlELP27d++1jB0Lzz1X6h45jlO1pBOVgw6CXXaJi0rHjtC3b3z9mjXbur+CdONgIGTE5Csqj2LCMl1ErhWRL0XYp4pgxx3hmmvm0acPfOMbMLuqbTPHcSJF1Z5Or7wyfbslS+D+++1zqhhIly5x91eHDjYo8tFH4YwzklsqAVu2FKTKcb6i8n3gMOBKoA54VkQWRtarCmHHHTfx1FPQrRuMHg3z55e6R47jVARBna4rrkjd5oMPzAoJRCWZpQJxUfn8cxMVEfj6181qCQL1RxxhN6rbb2+6bQHKVOUlKqp6v6q+r6q3q+qJQDfg9Gi7Vhn06gVPP23/y8MOg9dfL3WPHMcpe1atsve2bVO36d696fd0ovLxxxYrCY9Rad8+bqnsvruNpD/pJBOgv/wFxo4tSAZYNtlfHUWkVbo2qrpFVWeEttkxis4VkigLSvbta8LSooVVnS5wcoXjOJVOICpt2mS/TSpRqauDxYuhoQH23DO+vH17aGy0cSrhMS077QSnnWaB4Z49c+56xm5m0WY9cJmIbJ+xJSAipwBlP6Iw6oKS/fvDU0+ZVTtqFLz3XiS7dRynGgkm1wpbKkuXpt8mlVVRVwcbN9rnsKh06BD/HBaVApNNSnEjcD1wq4icKyL9EtvEyreMFZEHgJWqWtqpx0rEXnvZHDkffwyXXFLq3jiOU7Ykur8WLjRf+v/8T9N2YdFJZ6kEJFoqAd/6Vv59zZGsYiqq+jEwHtgRqBeRDbHaX8tEZC3wFjAWuFBVHyxcd8ufoUPLYkZPx3HKjM2bQ18SRSWoHPz97zfd6Hvfi3/ORlT22CP+OSjPMmFCfCriIpB1oD5msdwK9AH2xwLzFwCHAD1U9TRVfbcgvawwhgyxmT3D00c7jlO7rFsHw4bBo4/Ggu+JorIuRVWrcMwlG1Fp3Tr+OZhW+Ic/zL3DzSBrURERAV4FOgJXA0+o6n2qOktV1xaqg5XI4MGWhv7KK6XuieM45cAPfwgvvwzdu4em/AUTlU8+sSrDychVVMLst5/diAYOzK/TeZJLSnFb4DmgPzAYGBBeKSJXisjdInKZiHRNtoNaYcgQe3/55dL2w3Gc0jNtGtxy7afcc/D1XPTDA+Cf/7TAK5iojBtnKb7JCItKqkB9kHo8enR0nW4Gucz8OBZ4FtgXeAj4OGH92UAPYAswBfhGFB3MBRHpiA3I3Bt4UlV/Xew+gI1X2mknFxXHqXVWrbI5tB5r902GT59mC78emnqqvj59qZSWoVt0KkulbVuYNw92263Z/Y2CXCyVR7A5VF4HXsFEBtjqGusUi7vsCQzLtSMi0lZEmpvf2w8b7X848LVm7itvRMxacVFxnNrme9+z4QX7ydzkDTLV3goLSSpRARg0qGkKcQnJJVC/WlVvVNXnVPV6YK6I/FxE9gEuBJbH2i0Asq4FJiItYpN7LQL2S1jXS0Qmi8g5InKbiOyVoY/zVXUzcCDwp2z7UAgGD7aYSmNjKXvhOE6pmDIF/vpXG17Qcm2eWTvZikoZkXcvY2NR/gcYgZXCPza0LtE1lo6uwHRg1/DCmPXzMPCAqk7GkgOmiEhal52I9Aa+A1whImlqIBSWIUNsIOQbbxT2OAsW2Ij+BQsKexzHcbJk82ZWvLOKM880A+KnP23GvsJCUqBJtaImK1ERkYUi8isRaWJJqOpnqnqdql6qqi/m0wFV/VBVFydZdRgwEKiPtXsN2ASME5HDReSpJK9Oqvquqk4CXqaEk4XlGqzftAlOPhlezPEsPvaYVWf48Y9z285xnMKg555Hl3478NmKRm6/PbdKLNtQxZbKwcBi4CoReUtE/k9ERopIIf/Kg4B3VHVjaNki4FBVnaqqhyV5rQq1XQaUrArXgAGWMp6tqMydawVEb701t+PMmmXvDzwAL72U27aO40SP3HQjAFf+ZBODBjVzZxUoKlllf6nqh8CNwI2xYPrRWBzlVhF5CngAmJogAM2lO5BY7XElsEuqDUTkv4F9gLuBx1R1RYp2ZwFnAdTV1VFfX59XB1evXp122969hzJt2ibq6zPXxL/zzt5AP6ZOXUV9ffbq8NxzB3LAAWtYsKATF1zwGb/4ReIknZVPpvPsNB8/x9Hw8cetOS72ediQadTXt6f9u++yf0K7t08/nfbvvkv3DFPHvv7GG1vHbjz7/PNsaZbZUyRUNe8X0B74JnAn8F4z96XAYaHvfwSeTWjzN+Dh5hwn8TV06FDNl2nTpqVdf8opqnV12e3r8MNVQbVlS9W1a7PbZulS2+Z3v1O96ir7PGtWdttWEpnOs9N8/BxHw5gxaj9EUP30U1t4+unxZcHrlVeSLw+/RFRvvTX+ff36kv5tiQBzNck9tVn2lKquVZtbZQKWzhsly4DEFOPOQIZSntkRZen7VAwZYvPsLF+evt2mTTbVdL9+li02N0X2YSKB6+uAA+C737VpFZoVFHQcJ2/eftsyvrby17/CDTfAu0mqV3Xu3HQMSjJEqjdQHxCrRvx9EZkZKyi5TkTeFJF7scB6lEwD+iXEbQYQC9w3F4249H0ysg3Wz5ljc+n84Af2/YUXstv/rFl2Xe6zD2y/PVx8MTz+OMyYkXlbx3Gi5e674Vyuiy/47/+Gc86xueMTBaRz5/jc8+nuQRUYU8ml9tckzP3UGfg/4FTgKGyu+mnAiSJyj4h0ybUTKQL+M4EGYGSszQCgAzZav9kUw1IZPNjeM4lK4Mo+7jgbFDtzZnb7nz3bjhHUpDvvPNh5Z/jJT/LqruM4mXjnnaTzuqvCnXfCOdvfse02K1bAjgnzFnboEC9b/POfJz9WoqVSTaIiIucAH6nq4ap6uarepaqPqurTqvqAqk5W1ZOBi4CLRaRdth0QkZ2BS2NfJ4jIngAxn91YYJKInAdcBhylqtv+R/OgGJZK585WcXrevPTt6ustn71rVxg+3EQl00DbzZvNwjnggPiyDh3g0kttFsrp05vdfcdxEunXDw4/fJvFr7wCr70GO++4adttGhqsblPA1KkmGIGotG9vA9qWLGm6nQhst13T7xVANtMJdwGeU9XHQsv+IiJDE9uq6jLgR+QwPkRtnMovVVVU9VRVXRha95aqTlIbCzNJK3Dyr0zlWjZuhH/9C0aOtO/Dh1sMpqEh/X5ff90KnYZFBeDss6FHD7NWMgmT4zh5kMS//Le/2f1/p/brt22/YUNTUenf396Dchtt2tgc8olT+1ZrTEVVV6hqYp5qdyyQnqx9o6rOjqJzhaQY7i8wUfnPf5JazEA8nhKIyoEH2nsmF9js2BnePyFXsV07K7P97LPwzDN5dzsv/vd/U1fwdpyKp8ksW3G2bLF4yuGHQ8uN65Lf/MPur1at7D0QleB74udEUakQ8u3xD4AjRaRHlJ0pJsVwf4GJiiq8mmL4SH29XTtf/ap9HzTIrOFMwfpZsyy+Fzz0hDnzTKuUXExrZdkySxS48sriHK/S2bwZHnwQVq4sdU+ctCxfbj/QO++MzwOfwMyZ5ln49rex2kwdO27bKCwqQdA+EKmwi6ttQmWpGhKVRzFhmS4i14pI1gUky4ViWSqZgvVBPKVLLL2hZUubWyeTpTJrlrVLds21aQOXX25W+hNP5N31nLjlFvuNzJ8fnxnVSc7ChXDQQXDMMXDVVaXujZOW//zH3m+4IaWo/O1v5iEYOxabwXH77bdtlM5SCWeGhUWlxiyV72MpxFcCdcCzIrIw/SblRbEslb59oVOn5MH6xHhKwPDh8O9/p55hdO1aCwwmxlPCnHqqHfvii+H66+Hmm81Ef+ghE5pnn4U338zzj0pgyxb4858tlgPFd7tVCo2NcPXVlgL+xhtmZT70kMe+yprght/YuK2ofPQR2qIFe9z8Q8485E3TkmwslUBUklkqiZNyVaCo5DJJ11ZU9f7Yx9uB22MpwQdG1qsqIt3cKrNnm3AkE5XGRisu+eUvb7vdSy/Z9ZhOVFq3hl//GiZMsFTjVBx4IJxxBhx/fPLfQjZMm2YDv267Dc4/37LPTjghv31VK6++akI/dy4ce6wJ/f332//mP/+xWnFOGZNMVF57DVHlonW/YuNzN4J+bD/oZD+kcKA+LFSQXlTC6yqEbLK/OopIq3RtVHWLqs4IbbNjuva1xpAh5hZKjPMlxlMCgmB9qrhKMJI+MUifyPjx5rNfvtzS6197zYTq+efhySfht7+Fzz4zUenRA846y4Qu1yfnm26y38y3vmUC+fTTuW1fzWzaZC6uffeFxYvhnnssmaGuDo46ytpMiWTkldNs9tsPTjyx6bING+x98+ZtRaV9+60fW3/+SXx9tu6vM86w92CUNNRMTGU9cJmIJDlT2yIipwC7N6dTxaBYMRWwuMqaNfDWW02X19fbuvBDDNgAxn79UsdVZs+28S91dZmP3aGDtevbF/bc025uBx0Ehx1mA34XLDAX3HHHWSzygAPsGr/77uz+to8+sgrJEyfa7+Gww0zA3nknu+0rhRtvtDhVLoK7bh0cfLBNS3DssSbq48fHk4N697b/v4tKCbntNvNBgpmRiRd+4INOYqmsXxX63rdvvG3YUmnXbttlgVAce6xdUL16xdfVQkxFbYrg67GKxOeKyDY1vmLlW8aKyAPAykoYT1KsmAokL9eyYYMF0hNdXwHpBkHOmpXe9ZULIjBihMVcli61eKSquc3+/e/M2992mz2Nn3mmfR81yt6rzVr5/e/hF7+w2FG2fPe79j+8/Xa7V3Xrtm2bo482UV+RtJ62U3AmTYJx42zgVzLSiMqM6aGBjj16JBeVyy6z9z594svSjTepBVGBrTM5jgd2BOpFZEOs9tcyEVkLvIWNfr9QVR8sXHcrk4EDzY0aDtbPmZM8nhIwfLil6SbWovvgA0tfzOT6yocddoDvfMeC+F27WtmioDxRMlThT38yUfpSLP9vzz3t91VNorJqld1z2raFCy7Ibt6a224zAbrsMjjppNTtxoyxc/zYY6nbOAXk9tvtPTzaePVqM0s3brTAOyR1f01/0r5ru3bmigjahkXlzDPh00/jP5BMVEGgPpc56huBW4E+wP7A6cAFwCFAD1U9TVWTlON02ra1m23YUpk2za6Zr3wl+Tap4irBoMeoLJVkdO5sAxlnzTLRSMXzz1uQObBSwP6mQw+1DLBqyWp68UX7W/78Z7M2vvlNu0+k4tVXrbLBwQenLusUMHQodO/uLrBScMstMObqgwB474VQiZRf/tLM0r/8JaWl8skHm3h5jlkqsuuuVt4imaXSsmV2FYkDEi2VagzUB8TmjH8V6IjNF/+Eqt6nqrNUdW2hOlgoihlTAfOdh0UlVTwlYNAgc8cmxlVmzbLrbN99C9ZVwAZyHXqo1RL74IPkbW66ydKlv/WtpstHjbKxKqkGfFYac2LO3NGj4e9/txJNJ5+c3IpbvdrOR6dOcNddme8lLVpYwP7xx1MOg3AKxK23wnNvWWmUm66Ii8ry92JurdWr40KRYKk8cs8aZHOs3U47WdtklkoQlM9HVKC6LRWgLfAc0B8YDDRJghSRK0Xkbpdpz4oAACAASURBVBG5TES6RtjHglDMmApYXGXpUgtsB/GUQw5J3b5Vq+SDIGfNgr33bpJ4UhBE4LrrzKoPSvKH+fRTy2KaMMGSAcJUW1xl9mxLnOjSxSzI//kfeOQRS9kOo2ruw0WLTFB6ZFlv4uij7UHXi4AWj8ZG+79OPKMNW7p05dtfjYvKX+5oDcA/H97IsnfWxzfYFI+hPHLPGvr2iIlMly5NLZVw9lcgJtnW7aqVmEqMscCzwL7AQ8DHCevPBk4Cfo25yZwQQbB+3jy7mNevTx1PCQgGQQYPQFu2bFuZuJAMGGCDJ2+/3dx1Ye64w/oVdn0F9O5t9fGiEJXNm9PHdYrBnDlNY1jnn2/jcC6/vOl5ufFGG13985+nf2BIZNQou5e4C6x4vPKKDSIeMQJa7NKLPTvFReUbY826mPuvjdxwrQnFmlWNfPJB3FKZP2M1Xx0eslTWrzdrBbZ1f+VCjYnKI8Aq4HXgFUxkgK2usU6xuMuewLAoO1kNhMu1ZIqnBBx4oD0cvfiifV+0yMaVFCJIn4of/Qi+8AU499y49R8E6IcOtdHhyRg1yp68g/Fd+fD66+aOHjVqJC1bmjtwhx0siaBnT4t9vvde/vvPhg8+sGSJ/faLLxOxv79/fxOXJUsseH/hhXDkkfGEn2xp395SsadMyS4OVS2xqlxYv94sxLUROdqDQsPDh2NBrVBtoSFDLY7x3xds5OjDTFQ2f76WM06Oi8r2rGLEsJClAuaGgOaJyq9+BV//un2udlFR1dWqeqOqPqeq1wNzReTnIrIPcCGwPNZuAVBxtcAKTdeuVuTx5ZctnjJkyLbz9iQyfLi9B8H68PTBxaJdO/jjH+0G/9vf2rLZs+1JL5mVEjBqlHkE5jQjufySS+w3NWnSYi65xDKvTj3VRv8fcYSN+7j//sz7aQ5B/8OiAnbfuP9+ezgdP97iKHV1ZtXlcx8YM8YGR2aKQz35ZOqyP9XMlCnw/e9bHCQKZs4092SfPthFHq6JFHsSat9yI0O/ZG6C7WU13zoqXtZ+/z1W0r1LyFKBuDAlG5OSLTvvHM9Iy2f7MiCvMi0AqjpHRBZhLq/uwLGhdYmusbJDRMYAY3bfvXjjNAcPNmF4/3178s9EXZ1ZCUFcZfZsc9cWu6TH179uGU9XXmkDjv/0J4ujJA4+DhO4f55+Oi6OufDss/Dww5aIM3z4YkaO7LtNm5kzLcB90UW57z9b5syx33WyxIiBA+1cTJhgD6TPPht/aM2V8Oj6vVPMRrRqFZx2mgnZ7bfHrd9aIHiwuu02S3VvLjNm2HUpgqXxfv55fGVgkm/cuPWzqHLioR/Aw7bqJxeujLcLRCWZpZLPHCjBgMnx46s7+ysZqvpZbAKtS4E2GTcoI4odqAezTt58M7t4SkB4EOSsWTBsWGmus2uvtRvnmWdaEPqEE+yJORVdu9rfm09cRdWeSnfZJb1gHHmkudhSFd6MgjlzzM2WmIwQ8O1vwx/+YOckH/EM6NnT/rfp4io/+IEle+y5J9x7b+HdYAsX2n3twgsLe5xsCKz0F14wN3Bz+OADq/iw9f/Vpk08HgKm3mCiEb643n9/68e6Vp/GA/eBqHwce5ZubhZNu3Zm9Vx3XUVaKs3qsYi0EpGJIjIHC+I7aQiC9dnEUwKGD7cbyRtvmMujmK6vMLvsAj/7mYnE2rXpXV8Bo0bZE2GufvB777Wb+VVXxR/akjF6tAl0obKmVM06THR9JXL++VbmprmMGWM3z2Qp3E8/bSnc//VfFrN57734jTZqli61OnB77WUp1H/4w7YlhopJEFc8/ni7x96RZBr4XAgs/xEjYgsSLZVgkpsNG5qKSjiAtzJkqQRPV8FMfIlpwfnQrZs9xdWKqIhInYhcAbwL3AzsBFTGXJclJBCVffaxAHQ2BIMgr7/eXL2lEhWwsiNDhpgrKJtkgVGj4uX9s2XDBrtpDhqUfiQ6WCHOtm3hn//Mfv+5sHixlU/JJCpRMWaMCdmjjzZdvnq11R7s398yy44+2qpQ//3v0R5/1SqrU7bHHjYw8PzzLQbYogVMnhztsXLhlVfs4WHcOEtouP325mUEzphhKftbXZpt2sQLR0JcVIIR9a0txThsqbByZdxSCdxdgai0idBpU+2iIiL7icgdQANwDiYo/YBrCtC3qqNfPwsOBv7zbBg82J7W//pX+17MzK9EWrY0q+Dpp7NzFX/lK7ZNLi6wyZPNNXHNNZndfO3aWezm8cez338uBEH6Yp3zIUPMIkx0gV16qVUR+etf4xlwRxxhohJFuvXGjfCPf/Rit93MOjz6aEvM+P3v7fo75hg7dlSZV7kSTlA5+WQT+1weVBKZOdMyF7caFIkikGipBHW7AlFp166pqAS+0eAERWGpBFSrqIjIiSIyE3gBE5HTgV1V9Yex0iw1mOSYOy1amJ/68suz36ZVK/O1f/653XB69ixc/7KhU6fsrayOHc3SylZUVq60ZICvfc3m+86G0aPNx/7229m1z4XZs+1+kypwHjUiZq1MnRofmzR9urnWv/tdqy4dMH58dC6wCy6AP/xhDwYNMiG96y57AAo47zwb7Jpt5eqomTXLvEF9+5q10qGDBezzYeNGK0bcJP6VKCpBwb2gnlfPnvbjDdxfO+9sJ2TjRvuBBttH6f4KqPJA/WbgHuAoVb1TVb2oRB7ssEO8ckO2BD+AUrq+8mXUKPOHp6uVFfCrX1m73/wm+/2PHm3vhbBW5swx6yHX/1dzGDPGHnifecbeTz/dbvC/+EXTdkcfbfeye+9t3vHefdeskKOPXsJTT9kDTCIHH2zJCtddV5oxMkFVbhETlOOOs787nwSNl182nUgrKoFFsnq1CUe7dqZqgVnYrZv5RTdtMtdY4B4LLBV3f2VGVe9S1S9jo+WvEJE/isjQwnbNCQjiKqV0feXLqFF2I6qvT9+uocHcLRMnNp2zKBN77GFp11GLyubNJobFPueHHGI3zilTzKJ96y2ra5iYfdapk7nA7ruveS6wa66xm/WECe+mdGmKWAr8Sy8VLjkgFStXmisu/EA1caLFf/KpQBAE6VOKSrhE/Zo1ceHYay9b1qqVZXc9+aSZdLlYKkceyfKvfS23DlerqASo6jxV/S5wOfBVEblNRCZQYenEUPyCks1h1Ch7Okss3FgJHHCA/QYzucB+/GN7v+qq3PYvYqnFzzzTNNbaXBYutHtEsYL0AW3bmuvvrrssjfvcc1Onn48fbw/VqWYIzcQHH1jl5ZNPhp13Tn/yJk60MVLXXZffsfIliGuFRWXkSHMF5+MCmzkTdt3Vtt9KWFTCI5JXrIi7uIJA6KZN8MMf2uelS5uKSqaYymOP8XqwbbZUu6gEqOpKVf2dqp6M1QA7hAqLq5RinEq+dOpkQdkvfKHUPcmd1q0tSyudqPz735YmetFF9oPPldGjTQCefz7/fiaSaiR9MRgzxsrx9O69bdHKxHbNcYH97nd2z7zkksxtt9/exOfee+Nj/IpBYBmF/w/bbWeZgY8/nrqCdiqCQY9NCItK+H6wZIldWK1bW7ZCwBFHxC/UsPtrzRp7yonSX1orohJGVZ9Q1WOA3SLoj1OFjBplLowloSkrGhtt7M0jj1ggeKedcq+ZFXDIIfa7jtIFNmeOiXn//tHtM1vGjjWX5623Nh2cnUinTiao+bjAPv3U0tTHjzcXYjYE9d/+8pfcjpXIzTdvO/lcKmbNsgoSickhEyeaizKX5IH337dY+9bxKQGpRGXLFvPLtm7d1C0G8SmAEy2VqMeW1KKoBKhqQ+ZWTi0SlMI//3x74Bs40Fxi/fvb0/YLL9gTeb5GY8eOlr4cpajMnm1B61L8pnfaydw0Bx+cue348SbWiVMkZOIPf7CMwlyEfOBAE/AbbrAbej4sW2alZq64InPboIpEsgSVgQMtLTgXF1jSeAo0dVcluwgDy2PFCliwwD4HohK2VDZtyr2AZCaqPPvLcfJi8GBLB330UZsp8otftJHhN99s7oiPP7Ysp+YwerQVYwyPT8uXDRtg/vzSuL5yJR8X2OrVlhQxZowNMs2F886zh/fEAZrZMn++vf/jH5ljYIsXm6stVdbjySdb8kBwn8/EzJmmH9skgqSyVAIC0dhpJ1MzaGqphN1dUYtKLVsqjpOKFi1MTNautcrCDzwAV18Np5xiT42pZr/MhShTi+fNs4fOSsi22357S1TIxQV2443wySfxeHMujB1r99N8A/ZBdeXPPoMnnkjfNlNV7hNOsAf5cFHfdMycadZnoBFbyVZUwgSisn69xVGCfURtWbioOE5yWreO/iEuzJe+ZL/zKERl9mx7rwRLBcwFtnRpfI6QdKxfb1MYHHpoPFU9F1q2tNktp07Nr7Dj/PlWVaJLl8zxkBdeMMsi1eDTnXc2Qb3jjszuuPXrLUU8adHPVKISCESywHsQY1m2zN4D4XFLpfpERUT2EZEbSt0Pp7gEqcVPPtlk1te8mDPHph1oknZaxhx1lN18s3GB3XILLF9uk6/ly5ln2n02n3pg8+ZZ7bvjjoOHHooP7UjGrFkWN0mXTDVxosWUMo2Deukluy62CdJDalEJYi3JLJUgFTNxwKOLSnmIioi0FZFm5/aKyPbAoUCEdRKcSmH0aBsUl++4jYBg+uB8psIoBdtvb3PeZHKBbdpkCREHHJDbdMeJdO9u8+vcfHN6UUhkwwbLAhw82FxXa9emjs1s3Gip5pmqSIwZYzqQKWCfMkgPqUUlULNkohKuYxNu46JSWlERkRYiMglYBOyXsK6XiEwWkXNigyz3ymKX3wT+UYi+OuXPqFHmsWiOC2zVKrvxVYrrK+Bb3zJPTLpCi3fdZcHvH/2o+YJ53nkWF/nb37LfZuFCSyUfPNiy9Xr0SO0CmzfPRCiTqLRrZ+6/+++PT2eSjBkzTAfq6pKszMf9lTgbW6EsFc/+ypmuwHSgyZC32Jz3DwMPqOpk4Gpgioik/I+JyFHAP6mwQZhOdHTubO6N5ojKiy9aKmuliUomF9iWLVZbbdCg3Kpkp+KggywR6q67st8mCNIPGmT3yvHj4bHHTJwSCYL02cR9zj/fYirf+EbTubYCVFMMegxIJSqBlZDMUglUOVjn7q+tFDB0mhlV/RBAtn1sOgwYCNTH2r0mIpuAcSKyCrg4ye4+B84A2gMDROQiVb22QF13ypTRo+1JfPlyc9PkSilH0jeHjh3tpnrPPXZfXLPG3EvB+8cfmwV2113RuPVELNh/881mfWRzL50/34QvGGx54omW2vzQQ5YeHGbWLPv/ZVNhYdAg+7uPPdbcclOmNNWBhga7HiIVFbD89UQXmYtKyS2VVBwEvJNQCXkRcKiqTlXVw5K8jlHVccBZwDMuKLVJkFo8dWp+28+ebW6SfOeaLyWnnWbj8371Kxv1/vDDFkt44w2LUUycGG39uOHDTbRefTW79vPmWZZecN/df38bv5TM2glXJs6Go4+GP/3J/u+nnNI0trTNTI+JpKr9FdzQU2UK9OplKWjhfYRFJYpyLRUoKiW1VNLQHUg0ilcCkeTjiMhZmPhQV1dHfabUkRSsXr06722d7MnlPG/ZAjvuOIJbb/2U3r0X5nys558/kIEDV1Ff/1rO25aa9u1h6lShRQtNeTN+7rnky/O5llu0aAscyC23LGLcuKVp26rC3LkjGDFiBfX1/9m6fMSIL3D33b156KEZ7LCDpe2tWtWSN974Mgcf/Db19VnWc8ESss46a1duumk31q9/nwsueBMR+Pvfd6dt2x588snz1Ndv6x2XjRvZt39/Phw5kvcXLCAoZLBh0ybaAK+//TbLM5ybfdatYwfg8/XrebG+njZ33cWWdu3YlLBdzud582ZGxj5WzL1GVUv+wuIgh4W+/xF4NqHN34CHoz720KFDNV+mTZuW97ZO9uR6nk8+WbVLF9XGxtyOs3y5Kqj+9re5bVcN5HMtb9mi2r276kknZW67dKmd22uvbbr85Zdt+Q03xJf985+27Omnc+6Sbtmi+l//ZdtfdZUtGzpUdeTIHHZiGqjau7e933575m0OPdTaDhuWtlnO53nLlnh/ygxgria5p5arbbUMSEwx7gykfxzKgUoqfe/kxiGHmBvozTdz267Y0wdXOiLmUspm0GVQnmXw4KbLBw2ygpHhLLBZs2zfySYMy6ZP11xjVYwvv9xiNvPmpXF9ZdoZZOfGKlSgvlLy2kOUq6hMA/qJSLh/A4gF7qNAK6j0vZMbwRCCbCvhBsyZYy7sffeNvk/VyvDhNpVzphL04cyvMCI2ZmX6dKsKACYqAwdaFeZ8aNHCZrM88kibTqGxMU2QPhnnnGM5ypkC9WEKFaivQEouKgnCETATaABzJ4rIAKADkMdcbymP65ZKlRJU0GjIsW72woWw227bzrLopCa4WWeqkjx/vlUoSFbn7fjjzb/z97/b++zZzZ86u1Ur29/w4Za+nFNJmuuvt1SyYIxINqJSqNpfFUipBz/uDFwa+zpBRPYEcx4CY4FJInIecBlwlKrmMH43PW6pVC+9etlvO1dRaWiwbCQne4IyKplEZd68bV1fAQMGWOXgu++26ZNXrGi+qIA9HDzxhLnnunbNYweZsr/CFMr9VYGUVFRU9UNV/aWqiqqeqqoLQ+veUtVJqnpd7H1OKfvqVA4tW5qw5Coqixe7qORK27bmLkwnKkF5lnRl9k84wcrr3HOPfY9CVMBK2OQdI8vF/dW+vb27qJTe/VUq3P1V3fTunVtMZe1a+PBDF5V8GDHC4lEbNyZfHy7Pkorjj7f3q6+2+/OXvhR9P3PGRSUvalZU3P1V3fTpk5ulErR1Ucmd4cOttHwQjE8kVeZXmL59bT+rV1vWV1ncm3NxfwWBOI+p1K6oONVNnz5WRSPbaW8DUUmcitzJTBCsT5VaPG+eucl23z39fk44wd6jcn01m0BUslG4wFJRLz1Ys6Li7q/qpk8fc7kszXJk0+LF9u6WSu7ssovV6EoVV5k3D/baK/O9+fjjTXjGjo2+j3kRWB3ZCEUgKql8gDVEzYqKu7+qm9697T1bF9jixebl6NGjYF2qaoYPT26pqJqopAvSB9TVWZ2ygw6Kvn95EQw8zMbcdVHZSs2KilPdBG6sbIP1ixfbNhVYv68sGDEC3nvPZmEMs3y5VUhOF08pWzp3tvdsLopAVDZsKFx/KoSa/Qm5+6u6ydVSaWjweEpzSDUIMgjSZ2OplB133glXXZVdiQW3VLZSs6Li7q/qpkMHG/CWi/vL4yn5M2SIBeMTXWBBRlhFWirdu2c/TaaLylZqVlSc6qd37+xEZd06c9O4qORP69aWCpzMUtl116bTlFQl7v7aiouKU7X06ZNdTCVo46LSPIYPt+mY16+PL8s2SF/xFNpSyavOTGmoWVHxmEr1EwyAzJQR6mNUomHECNi0CV56yb4H5Vkq0vWVK4UUlYYGWLQo+v0WiJoVFY+pVD99+th0t598kr6dj1GJhsRBkEF5FrdUmknv3hXlP6xZUXGqn2wzwBYvtoF5PXsWvEtVTV2dzWUTxFWyKc9SNXigfisuKk7Vku28KosXmwB52abmEwyCDAY9ZlOepSpwUdmKi4pTtWQ7ANLHqETH8OGWSdfQkH15lqqgXTt7v+CC0vajDKiFf7dTo3TpYg+Q2VgqRxxRlC5VPcFc8DNmmKiUTR2vQtOiBWzZUpFzykdNzVoqnv1V/YhkHquyYYMVnfQgfTTsvbcNPH3gASvPUhNB+gAXFKCGRcWzv2qDTPOq+BiVaGnZ0mZafPBB+14TQXqnCTUrKk5tkGkApI9RiZ7hwy2VGGrMUnEAFxWnyunTBz76yKYLToaPUYmeIK5SE+VZnG1wUXGqmmCsSiprZfFiSyXu1atoXap6DjzQ3t31VZu4qDhVTaaxKg0N9kRdE2mvRaJLFzj7bDj55FL3xCkF/lNyqppMohJMzuVEy+TJpe6BUyrcUnGqmp49zb2Vzv3l8RTHiY6aFRUfp1IbtGwJu+yS3FLZuNGmv3VRcZzoqFlR8XEqtUOqAZDvv281qlxUHCc6alZUnNoh1QDIIJ3YYyqOEx0uKk7V06ePubmCAXkBPkbFcaLHRcWpevr0gc2brcZXmMWLrQ7gLruUpFuOU5VUraiISAsR6VDqfjilJ9VkXQ0NJiitWhW/T45TrVSdqIjIxSLyJvAisKXU/XFKT6qxKj5GxXGip6pERURaA12AvVR1H1VdV+o+OaUnlaXiY1QcJ3rKRlREpK2INDe/dzdgEPCuiHwzgm45VUD79tCtW9MBkJs2WUqxi4rjREvJRSUW+5gELAL2S1jXS0Qmi8g5InKbiOyVbl+qulBVjwS+DPw6Zrk4zjZjVZYssYn6XFQcJ1rKofZXV2A6sGt4oYgI8DBwmapOFZHpwKMisgdwKHBxkn0dq6qrVHWRiDwGdAI+Lmz3nUqgTx947bX4dx+j4jiFoeSioqofAsi2U3EeBgwE6mPtXhORTcA4Vb0PmJq4gYi0CX1dr6ouKA5g4vHPf9oIehEfo+I4haLk7q80HAS8o6obQ8sWYVZKKi4RkQdF5CTgwYL2zqko+vSBdetgxQr7vnixicuuu6bdzHGcHCm5pZKG7kBitceVQMqhaqr682x2LCJnAWcB1NXVUV9fn1cHV69enfe2TvZEcZ5XreoK7MX998/li19czaxZX6Rr1x2ZMeOFSPpY6fi1XBxq4TyXs6g0ApsSlkViWanqTcBNAMOGDdORI0fmtZ/6+nry3dbJnijOc6dO8JOfQLduwxg5En72M+jfH///xfBruTjUwnkuZ/fXMiAxxbgzsDRJ25zx0ve1ReIASB+j4jiFoZxFZRrQT0TCfRxALHDfXLz0fW2x007QoYOJSmMjvPeei4rjFIKyEJUE4QiYCTQAI2NtBgAdgCkRHdMtlRpCxMaqvPuuFZbcvNlFxXEKQclFRUR2Bi6NfZ0gInsCqKoCY4FJInIecBlwlKquieK4bqnUHsG8Kj5GxXEKR8kD9bFxKr+MvRLXvQVMin29LsrjisgYYMzuu+8e5W6dMqZPH5gzx8eoOE4hKbmlUircUqk9+vSxcSoLFtj3oNCk4zjRUbOi4tQegYg89xz07Alt2qRv7zhO7tSsqHigvvYIYihz5ng8xXEKRc2Kiru/ao9ASBobPZ7iOIWiZkXFqT169oTttrPPLiqOUxhcVJyaYbvtbE56cFFxnEJRs6LiMZXaJHCBeUzFcQpDzYqKx1Rqk0BM3FJxnMJQs6Li1CZ77AGtW/sYFccpFC4qTk1x4YUwYwa0a1fqnjhOdVKzouIxldqkUycYOrTUvXCc6qVmRcVjKo7jONFTs6LiOI7jRI+LiuM4jhMZLiqO4zhOZNSsqHig3nEcJ3pqVlQ8UO84jhM9NSsqjuM4TvS4qDiO4ziRIapa6j6UFBH5CFgJJAZXdshiWVfg48L1LmN/CrV9Nm3TtUm1LptzmqxNpZznXLfN1L5Q5zjZsko5x7luX8prOdmySjnP2WzbR1W7bbNUVWv+BdyUzzJgbin7WKjts2mbrk2qdVme02RtKuI857ptpvaFOscpzntFnONcty/ltVzJ57k527r7y5jSjGXFornHzmX7bNqma5NqXTbntJTnuLnHz3XbTO0LdY6zOXYhqZVrOdvjF4piXstbqXn3V3MQkbmqOqzU/ah2/DwXHj/HxaEWzrNbKs3jplJ3oEbw81x4/BwXh6o/z26pOI7jOJHhlkoREJG9S90Hx8kXv36dXHBRKTAisj/wQqn7UU2ISCsRuUJEjhGRH5S6P9WMX7+FRUS6i8g/RKRBRH5W6v5EgYtKgVHV2cBHpe5HuSMibUUk25o5ZwBvqeoDQEcR+XIBu1bT+PWbOzley4cA44G9ge+ISOfC9aw4uKikIMcLw8kTEWkhIpOARcB+Cet6ichkETlHRG4Tkb1iqw4E5sc+vwyMLl6PKxu/rgtHntfy/araqKqrgNeAdUXuduS4qCSQ54Xh5E9XYDqwa3ihiAjwMPCAqk4GrgamiEhLoDuwOtb0c2DbUb1OE/y6Lgo5X8uqujHWphvwlKpuKHKfI6dlqTtQhmS6MC5T1akiMh14VET2AA4FLk6yr2NjTyBOClT1QwA7vU04DBgI1MfavSYim4BxwAqgY6xdR4pX9qKSyfm6VtXGEvSzYsnzWr4v9j8YA/y6aJ0tIC4qCeRzYajqfcDUInazFjgIeCd4kouxCBPwaZgPeh4wCHi6+N2rLPK94RWxi9VMumv5PuAY4B5V3Swiu6rqe6XoZFS4+yt70l0YKRGRfYFuInJ4ITtXhXRn24J2K4FdgJuBPUVkPKCq+kyxO1dFpL2u/fqNhJTXsoicA/wOmCUii4ABxe5c1Lilkj3pbnIpUdWXgA6F6lQV0whsSljWAiDmlvlR0XtUnaS9rv36jYR01/JkYHLRe1RA3FLJnpQXhlMQlmHlt8N0BpaWoC/VjF/XhaemrmW/eLKnpi6MMmAa0E9EwtfoAGK+fycy/LouPDV1LbuoZE9NXRjFJOGcBswEGoCRsTYDMDdMqUvjVxt+XUeIX8seU0lKFhfGM9V+YRQLEdkZGyEPMEFElqjqQlVVERkL/ERE9gT2B45S1TUl62yF49d1YfFr2fAqxQmELoxfALcAv1HVhbF1uwE/AWZjF8YfVXVOibrqOFnj17VTLFxUHMdxnMjwmIrjOI4TGS4qjuM4TmS4qDiO4ziR4aLiOI7jRIaLiuM4jhMZLiqO4zhOZLioOI7jOJHhouI4juNEhouK4ziOExkuKo5TZESklYhsX8TjdS3WsRzHRcVxioiIdMMmGNsoIt8QkVdFREVkYJK27UTkExFZIyKnNeOwPUTk3GZs7zhZ46LiOEVCRFoBNwG/U9UNqvoo8BA2Udb5STY5Prbu36r613yPq6qvAO+KyCn57sNxssVFxXGKxyXA06oax7BExQAAArNJREFUnr53E/AAMFFEOiW0Pwp4AhOWZqGqj2Dl2Hdu7r4cJx0uKo5TBGJzmXwHeDzJ6uuB9sApofZDgPnA5oT9nCYic0XkFBGZJyIfiMiZ4e1E5BcicqmIPCUi/UKbzwYmRfdXOc62uKg4TkSIyHEislhEviQij4nIzaHVg4Duqvpmkk0XA48C54mIxJadDvw5SdspwFBs2orBwNXAZBHZXUS6A5OBn6rq1cAKTMgCFgDHNeNPdJyM+MyPjhMds4ANwAjgLKBHaF1v4LNkG8X4I+bqOlxE/gV0VNWlcY0xVPWj2LJpsUV/AC4HvgZsD8xS1cBdNgkIT5i0Ctgt9z/LcbLHRcVxouMQYDvgMVVdArwfWteGBFdWAk8C/8EC9l8A7sjmgKraKCLvAO0wwVgbWrc+ofl6THgcp2C4+8txouPrxAUlkfeAHVJtqDYF6/WxfYwDnsnhuB2BhcCH2FzzWxGR3UNfOwHJ+uY4keGi4jgRICLbAYeTPA4C8AqwKUn2VQdMFMDmjl8LPKHxeb5bA62S7K9D7Lj9seywqcB9wGAR+b2I9BeRbwH9Q9t0B17K5e9ynFxxUXGcaDgQ+EhV5ydbqaprgPuBg4JlIjIWOBG4TER2V9VVwI2YuCAi38QsjyEicmrCLs8WkSuAnwLfVNXNqjoPC/AfAzwP9FbVx0LbDA/27TiFQuIPRI7jFBIR6QP8n6qObeZ+FPiCqi7OYZuewO+B8eo/eqeAuKXiOEVCVRuAW5tZciVAMjeJNbSR/N8DznVBcQqNi4rjFBFV/QewQES65LO9iJwd+3hKbFxKNuwLXKGqH+VzTMfJBXd/OY7jOJHhlorjOI4TGS4qjuM4TmS4qDiO4ziR4aLiOI7jRIaLiuM4jhMZLiqO4zhOZLioOI7jOJHhouI4juNExv8DBA3X6kXndSYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "frac diff for r = 0.1-100:  0.002526427174667445\n",
      "frac diff for r = 50-250:  0.15306941576688715\n"
     ]
    }
   ],
   "source": [
    "# Plot relative difference\n",
    "frac_diff1 = []\n",
    "frac_diff2 = []\n",
    "abs_diff1 = []    \n",
    "abs_diff2 = []\n",
    "\n",
    "for i in range(nr1):\n",
    "    frac_diff1.append(np.abs(xi1[n][i]/xi_toolkit1[n][i] - 1.))\n",
    "    abs_diff1.append(np.abs(r1[i]*r1[i]*(xi1[n][i]-xi_toolkit1[n][i])))\n",
    "    #print r1[i], xi_toolkit1[n][i], xi1[n][i], r1[i]*r1[i]*(xi1[n][i]-xi_toolkit1[n][i])\n",
    "    \n",
    "for i in range(nr2):\n",
    "    frac_diff2.append(np.abs(1-xi2[n][i]/xi_toolkit2[n][i]))\n",
    "    abs_diff2.append(np.abs(r2[i]*r2[i]*(xi2[n][i]-xi_toolkit2[n][i])))\n",
    "\n",
    "# esthetic definitions for the plots\n",
    "matplotlib.rcParams['mathtext.fontset'] = 'stix'\n",
    "matplotlib.rcParams['font.family'] = 'STIXGeneral'\n",
    "matplotlib.rcParams['font.size'] = 14\n",
    "\n",
    "plt.plot(r1, frac_diff1, 'b-')\n",
    "plt.plot(r2, frac_diff2, 'r-')\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.xlabel(r'$r$ (Mpc)')\n",
    "plt.ylabel(r'$\\Delta \\xi(r) / \\xi(r)$')\n",
    "plt.grid(which='major')\n",
    "plt.title('Relative difference') \n",
    "plt.savefig('benchmark_rel.pdf',bbox_inches = 'tight')\n",
    "plt.show()\n",
    "#print frac_diff1\n",
    "print(\"frac diff for r = 0.1-100: \", np.amax(frac_diff1))\n",
    "print(\"frac diff for r = 50-250: \", np.amax(frac_diff2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEoCAYAAAB7ONeTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2deZwcZbX3v2cm+zqBBLJvJAFCIIkkXLZ4oyyKgiwK6EVlkyDwul1F5dUXVC5e5YrovSCIEBC9Liib7EFgAEnYQwCJYCBknZCE7NvMJPO8f5x+0jWV6u7q7upt5nw/n/l0d3UtT9dUPb865zznPOKcwzAMwzCSoq7SDTAMwzA6FiYshmEYRqKYsBiGYRiJYsJiGIZhJIoJi2EYhpEoJiyGYRhGopiwGDWDiFwiIqdk+O5MEblbRH5aguP2E5GvicgcEfm/Se8/j3YcLiI3iMjcwLJPiMj7IjI0tO50EfmeiDwuIi+KyJDU8q+IyGUislxErin3bzA6ByYsRi1xAXBRhu/uASYAvUpw3M3ALcDhQLd8NhSRngm24zlAgKCIrACeAHYEjlkP/Bm4FvgEsAjoISLHATOdc/8JfJ48f4thxMWExagJROQwoB9wnIiMC3/vnGsG1pTi2E7ZBKwrYPNfJNkOYFVo2UvOuU8554Jt+yAw0jm30Tm3xTn3aefcYuCzwMbUdo87576UVNsMI4gJi1ErnAOcDLQCX6xsU+IhIl9E212OYwXv5WEZVsu03DASxYTFqHpEpC/QzTn3GnAncK6I9Miwep2IXC0iG0TkPRG5LLCfBhH5LxE5V0RuFpH5ge+6ishVIvITEblORJ4VkZOytOlQEVkkIo2pz2NF5EYRcSIyWkSGoW4oUvs8L/V+QOrzrSLyqojcks1dJiLDROR3InKNiFwFTAt8N0REvisi/wRGppZdBJwVOO5PUuv9BNgfmJZadk629ohIHxE5W0SeFpFDUufjdRGpF5GRIvLfIvIbEXkjdb7rRGSQiFwuIm+JyBQReVBEtojIHSn3nG/3mNQx/yMVA7pcRKSQ82NUKc45+7O/qv4DLkRjAwBHAw74fMR6jcAy4ExgEjA7te6/pb7/PvDlwPo/DLz/DfDTwOePA23ACYFl7wLfC3y+HWgMfP5Q6nijU5/PIeXBCqzzJ2Bo6v0QYCdwdYbf3RNYCHwisOxB4N3U+/7AxcFjZjpu4PzcFqc9wN7A11P7/gFwLHATUJ9qQ5/UNtNT61wMdEFFzQE/BgYCR6bO4ymp9YcBLwANqc8fSa1/fL7nx/6q96/LHkpjGNXHDOfcLwGcc38TkVfRIP7tEevOc879EUBELgSOA74G/A7oDswSkbucc8uBn6fWG4fGH47wO3HOPSAiLwPfAx7K0K620OesFV1TcaIjgC+nHtBBO/veGTY5FxgA3BdY9jwwMdXGjSLyRrZjFtoe59z7IvJSatkdzrnXgb+KyJnAKOC7gW3+BuztnNspIitSy25wzq0F1orIKuBAdIDFN4GHnXMbUuvNAc4Gni3g/BhVigmLUdWIyFRgsojcE1wMHC4iU51z80ObtPg3zrlWEXkM+FRq0c+A04GFInIFKWEBDk29bg3taz7a6SXFB4Clzrlvx1z/w8Ayl3p8LwFx27Mj8H4qsCCP3wD6P+meev9B4Jf+i9Rvux1ARPI9P0aVYjEWo9o5D/iwc+4U/wccgwbxMw09DrKGVMfonFuFxih+B/wX8JiIdAd2pdYdHtp2LeqKSYrewOjwQhHpEgq+B9dvSPD4xbbHbzMmYpu4Q5e7kooHJdQeowqxf5ZRtaSC9oOdc+2GEac+PwT8m4j0y7GbYcDjqf0d55xb75y7EPhXNC/lBNS91IbGb4IMBeaSGUf7eyjcubrUcb1f55/AkIhBAV8m/UQfZBEwRkRGZGlDMeTbHr/NYSIyObT86zGP+Q/g8yKyO99IRPqKyLEFtseoQkxYjGrmAmBehu8eQJ9wLwgsawN2jxYTkeGo6+W7qUWfFJGRoLEa4E1gpXNuKXAzcKGI7JXatj8aWL4isP+uaPDasxiYIiIHpjr/01PLR6defW7JASIyBXgktc1tInK+iMwQrRSw1Tm3PeI33oS6/W4Uzf7vnfo9e4vIYamn+K6pdYOi1i31G8KdcTfad9AP5WiP7x+CLvPfAluAv4hWO5gpIr8GFoTWlcA23QL7+ikq9k+LyFki8ingRuCZGO0xaoVKjx6wP/uL+gM+B2xCO5spoe8ORDPLHbAB+Fxq+QzgsdR3/4OOCjsosN1t6BPzV4HvAF8LfNcFuBLNYr8KzbT/UOq7PugAgDbgJeCDqeV7o8HlbWhg+oPA0+gIqd6pv5eA5cA5qW0mocHuHahFclGO83Ai8Bqa/X93qo2/R0etHZR679C4xfhUG55KLbsqdbxe6EixbUATOnKrf7b2APuhQ7v9voPn8YPAK6ltXgc+mVq+D/C/qW2uRt14s1KfX/D7SP1vF6MCdS8wIrDvvM6P/VXnn6T+mYZhGIaRCOYKMwzDMBLFhMUwDMNIFBMWwzAMI1FMWAzDMIxEMWExDMMwEsVKugADBw50o0ePLmjbrVu30ru3lTIqJXaOy4Od5/LQUc7zSy+9tNY5NyjqOxMWYPTo0bz44osFbdvY2MjMmTOTbZDRDjvH5cHOc3noKOdZRJZk+s5cYYZhGEaimLAYhmEYiWLCYhiGYSSKCYthGIaRKCYshmEYRqJ0amERkZNE5KaNGzdWuimGYRgdhk4tLM65+5xzs/r371/pphiGYXQYOrWwGIZh5MXWrXDCCfD225VuSVVjwmIYhhGXhx6Chx+Gb36z+H3t3AmXXw6bNhW/ryrDhMUwDCMbTU1w9tmwYwfUpbrMtraCdrVrV+DD//4vXHkl/L//V3wbqwwTFiNRVq/WBzEjN0uXwuc/D2vXVrolnZytW9PvN28O9f7A178Ot98Od94J9fW6LLxOmMsvh/PP3/1x50649VYYPx4WL+6lC7ds0deWliJ/QPVhwmIkxpYtMG4c/PKXlW5J9bN9O5x6KvzmN/Doo5VuTSdm4UIYOxb+8Aft/fv1g0suab+OiL46lxaWXBbLlVfC7Nm0tcGf/wwHHwznnQd77w0tLalut7VVX7t1S+73VAkmLEZiLFigD3wvvVTpllQ3zsHFF8PLL2uf9cYblW5RJ2bvvWHgQPjiF9Om9uzZ7dfx7i/n0u9zWSwppk+H00/Xze66C55/HvbfP2WpeGHp2rXIH1EAzukTTczfkS8mLEXw4IOwdm3lnzb+9Cc45BB1AVeS+fP19Z//rGw7qp0bb4TbblPX+oQJJiwVZZ994NOfho0boblZl4WtkWBcJUaMZe7c9Pt16+DXv4ZXX4VTJ/0TebIx/WUlheX+++H44+Haa0uyexOWAtmwAc46C77xjcmsWVPZtvztb/DaazpYpZK88oq+vvVWZdtRzcydC1/5io5YveIKmDjRhCUu992n52v79oR33K+fvvpE6fBTvHeFtbXljLEsWwbBivhvvqlxtPp69CniQx9Kf1lJYVm2TF8XLSrJ7k1YCqShAe65B5qaenD88bB+feXasnKlvt5xR+XaAGmLZfXq9D1qpGlqgk99CkaM0AFB9fXaUS5a1CHjt4nz9NMaEnn99eT2OXs2/G1BX/2Q6SbOwxX2yCOwszVtzWQNn1RSWEqMCUsR/Ou/wpVXvs4bb+gT6ObNlWnHihX6+pe/wLZtlWlDa6ve8OPG6Wdzh7WnpUV97Rs3wt13w4ABunziRHXtl+jBsUOxdKm+LliQ3D6vvhr+59cpi2XDhuiVghaLF5QMrrBHH4UDBmd/qpr5oQ/BmjWVEZbXXlPXgv9NJcKEpUgOO2w9d9wBL74IJ55YmY595UoYOVJHTT70UPmPD/ok2dICZ56pnzuiO+wXv9Cn5kL4+tfhmWfg5ps1HuaZOFFfzR2WG++9SUpY2trg3XdhQ5taLLvW5rBYcghLWxs89hicfPh7uQ8+d25lhOWQQ2Dq1JIfxoQlAU4+WYeNPv20DiH1McBy0NamwnLGGRqHrJQ7zLvBPvlJfRjqaBaLc/CNb8DXvpb/trffDtddp9t+5jPtv9t/fxsZFpekLZZVq/RenXCoWiyP/imDxRIUFj9yLMIVNn8+vP8+/Oc9B+Y+eGurucKM3HzmM/o0OmeOPrX7a6bUrF2rxxo5Ujv1++9vn+9VLl55BXr21AeiUaM6nsXy/vsaNH7pJR3hE5fly+HCCzWge/XVe37fsyeMGWPCkoudO9OxxAULVOiLZfFifT3jCyosj98VsFi2bdPIO0S7wiKEJa98pKCwlItf/apsh+rUwpJ02fzzztMn03vv1ZEgJRoi3g5/sw0dqoK2bRs88EDpjxtm/nwVlfp6zS7uaMLi3TCgGdRxufFGfSqePRu6dIlex0aG5WblSu3XDz1US2u9+27x+/TCMmSCusL2loCwnH46HHCAKpq3WFpacgrLwQcHFvTsmX4fXj8oLOXoKDZvhlmz0p+95ZWEQkfQqYWlFGXzL7kEfvxjTeTNpwMqFB+4HzYMjj4aBg+GP/6x9McN4pxaLFOm6OcJE9QVVqJrtiIsX66vY8eq2zPOKK7mZn1IPPFEtUoyMXGiPhxbKZzMeGE/6SR9TcId5oVl+ES1WD52VMAV9thj+trSEi0soRjLtm067P+4YwMXff/+6ZsgXGjy7LN1FAeUR1jCSlzipLdOLSyl4tJLoVev8jyFeotl2DC1Fk4/XRM3yzlC7d13dbSTjwmOH6+fK53fkyS+Y/vOd9Qtdt99ubf585916PX/+T/Z15s4Ufusd94pvp0dFR9f+djH1DOVlLAMGQI9BvYBYNLQtMXSVpfKV2ltTbvCooTFOWhu5umn9evjjwmIxKpVaqbOnRs94swPby6VsDz/fPpEmbDUPiIwfHj6KbeUrFihxxs8WD+fcYZeM/ffX/pje3zg3gvLhAn62pHcYcuWaR/xuc+piMexRq+/XkX22GOzr1ftI8OqwfL0wr7//npOkxKWMWPQf2yvXsiGtLDsaE0JS0tL+gQ0N+8ZvL/tNujRgxf//C7dusGMw0Nxk7Y2HRKYzd1eKmH5l39RN0Jz857C4gOxJTKTTVhKRLmEZeVKHQ3mB5YceaR2fOV0h73yinoLJk3Sz15YOtLIsGXL9Lx27arxs4ceSluLUcyfD/PmaU2wuhx32QEH6Gs1CsumTXDQQfD971e2HUuXqmepXz+YPDlhYQHo3budmd+yM+D+8p1vVIzlL38B4P1HX+aoo6BXt4iOevNmrYyciVK7wh5+eE9h8dZSiSwXE5YSMXx4Ov5RSlas0MC9p65O3WEPPVS++YPmz9fOsVeqGvioUfoQ2NEslhEj9P255+qD6G9+k3n966/X83HOObn33bevjuqrRmH58pc1R8mHAyrFsmV6jkCF5Z13iru+W1t1n7uFpUuX9oGzVOmWjWta0kH2KFfYXnsBsGvJMo2vRI308qM3MlFqYdm8Gd4L5dasW6evJiy1xbBh+kSbzzXT2qrB3nys0xUr9FhBzjxT74F7742/n2IIBu5B79H99utYFsvy5fqwAOqKmTFD+4ooN9G6dVqy5bOf1dI/cajGkWF/+pMWUBwxQodYZ0pMLwdLl6aFffJkfc1n2HeYZctUG9oJSyABrXd/HcJ37dWB0VvNzRmF5ed8lX/b8IvoUvjNzdk78CSFxbk9S9O0tOx5fBOW2mT4cBWI1avjb/PIIzoi8PHH42+zcmV7iwXUtTpyZHmSJdeu1U43nMw7YULpLJbt27VjmD+/PKVQ2tr0N/qODdRqeeut9pVsPbfeqvdreFqPbEycqJZBOQYIxWHFCs2/mT4dbrlF+6uo31oughaLf4gpxh3mR4TtFpb6+nbC0rW7Wiz33NHClg1ZXGG+Ng8w8vUH00+FYWHJljWd5D/94YdV7H796/SyKGExV1ht4p9u84mzeDdo3NFBzc068ipssYhoEP+RR0pfHDMcuPeMH68WS4EzuO5m0SL41rcOZupU7Vh69dK/kSPhAx/QvIFSJ4SuWaP3ZlBYTj9d3fLhIH5bm5Z+mTGjfemWXEycqPf4kiXJtLkY2trUhdfcDL/9LRx1lMaWnnqqMu3Ztk1H4vnzP2yY9p2JCksGV1jP+hYWvhpwhYWC9653n92byPBh0RZLS0v5hMUPnwv6YFtbzWLpKBQiLP6aiCssq1bpa1hYQIWltTWzO2zxYvjEJ7SDPuIIrbr71a/Cf/0X/P732on4mVOz4UvlB11hoBbLjh3Fx5luuw1efHEvhg+HD39Yg+E//CHcdJMO696xo/QuNz8iKSgsffroOf7jH9sL28MP6/8vH2sFqmtk2H//N/z1rzpVx4QJKuTTphVeJ61Y/Pn3FotI8QH8xYtVO/x9GnaFeWH59GktrFyS2RX2XlPgyWnIkLSwdO+eXr59uwpSplLHUcKyfHlhvsdgUqbHXGEdB3/B5tOx+qdV/zSVC7/vsCsMtCMYM2ZPd1hrqyZwHnQQPPGEJlX27g1//7u6PL75Tfi3f9PKzb6gZDbmz9cOd++92y8fP15fi3WHzZkDEydu4r77VGR+8hO47DK44AKdDyeJY+QiSlhAKy1s2aL5Kp7rr9eh36eemt8xDkyVl1q4MPe6u3aVbmDI66/Dt7+tiYgXXJBePmMGvPBCCeZCiYF/4Aqe/8mTtVBvoQ/7ixerUO2uhpBBWD7/6Ra61UW4wlLC8tY/Ag2oq0tbNEFh8fv1o1vCRP2IESPSN1E+RA0eiLKY/Ag4E5baYuBAdR8UYrHkKyxRFot3hz36qLoRAJ57TgXn29+Gj35UO7Hf/U6fThcu1Gtt40Z9ar7kEh1Zlqst4cC9J4lclvff16rR06ati/zel+j3JZ1Khf8f7n66TXHUUXrv+wE/b7+t5+zCC/OfxryhQR9441gs11yj5zfpUX8tLcJZZ+mw3ptvbl9ZfcYM7bOefz7ZY8YhbLGACsv27YVbq+2GGoMKSYQrbEDvFvYbpZ315vf3jLG8/WZAFJqbs89jH2VNBPbFrl3tb7i1a+P+nDRRJSGiXGEeE5baoq5OO/x8hMVbLHFdYcE6YVGccYY+QN1+u2Z/H3GEdtb33KPzb4c7StA8gQMPhG99SzuWm2/OfHxfpy+qCvfQofqAVoyb6rHHNGg8fXp0oKh3b32wK4fF0q0bDBrUfrmIBvGfekpjQTfcoP1RsCRTPsQdGfb73+u5T3JeEoDZs8fw6qsqlPvs0/67o47S31uJOMvSpXrs4ANUsQH8PYQlg8VCSwtjhqlYLP5HczsRaGmBd98OCEtLS7QrzJPJYpk9W5/Q/uM/tGbQ228X9qMgs8ViwtJxyCdJsrlZZxjs108D7nHqYq5YoR1e2A3lmTpVn+r//d81oPylL6llcvLJufc9YoSWz5g9O3MR1ldfVY9AlLDU1emxi+n058zRJ/n9989cn2bChNJbLMuW6f8yKtHx85/X5b/4hZ6rU0/NLPS58MKSLdP9nXfScS3/mgSPPw533DGCiy6Cj398z+8HDNCBEpWIsyxbpu7FoBFw4IGqBYUIy7ZtmtaRVVgC9cG6ojfAhjUtrFyacnW1tfHss7CrNSUsffu2z8yPIyxduugNv3Wr3kS33KLLiwm0RVksJiwdi3yExa83Y4a+xnGH+aHGmSaDE9GYycyZ6gb7+c/1+o/LrFk6QCBTeZhMgXuPL0ZZCM6psBxzDNTXZ+5p999fxauUZUeCyZFhhg2Dj3wEfvYzfSDIVRcsGxMnaswm2zXjExV79UqPyEuCb30Lhg7dzk9+knmdGTO0mkC5i2UGc1g83buruBQiLH705R6usOAQRh98CYwE613fTONjaYvl0Uehi6Q+9+zZ3mKJ4wprbVWz2+N9fvnkKISJmuOltbW9aN5+O/z0p9o5/OxnhR8rCyYsJcRn38fp9Hx8ZeZMfY3jDotKjgxzwQUapJ8+Pff+wpxwgv6Gm26K/n7+fLUoRo2K/n7CBP0dhUw78eabep8df3z29SZMKH3By3AOS5jzztP/8aRJ6QeDQogzMuzuu1XIjz46OYtl82Z4+WU45pjVGb01oL9ty5ZkLaU4BHNYghQ6MmyPocaw55wG/mktUN5++D7NvP3PXbuXP/oojBia+tyjR/sYS1xXmHe5BSlGWLzFEjx+2GIZMkRnnfvxj7X4XQkwYSkhw4bp/3NddOy5HT6+4oUlrsWSS1iKoUsXOP98zYeJmv/CB+4zWUzjx+vDXiFzZ8yZo6/HHZd9PT9IoFTuMD8CKyoe5TnpJI1BXHFFcVOJ5xKWpiZNUjztNPWcvP56vPL9uXjuOX1YnzQpu//Vi2axcZYtW3TYeJz7wrloiwVUWFasyD/GHUtYvPUSsEIG9d5Gv16p/JUdO3jhBRg3epeKQ/fuuYUlKngfJSzh8iv50NKi+wxaLF5YjjlG/dOHHVb4/mNiwlJC8sll8RbLwQfrqJxcwuLcnnXCSsF552ln6d2/np07NcaSbfrsYopRzpmjwpRtHhNQVxiULoD/3nv6W7NZLN2761wcn/pUcccaOFAHCGQSlnvv1f/7aaepoLe2xhuenItnntH/8YEHZh9mNnSoluopNs4yd64OdIgz9cC6dTr6K5PFAvlbLYsXax+/776BhdmEJeUKq9u6hSP/RYVFdu1C2nYyZlRKWLp1a59AGeUKi7JYojKIs1U3zUVrq4pKULB27NDlM2bozdivX+H7j0mHEhYRGSwid4nIEhGpcD3W/IRlyRK1ULt314EhuVxhmzZpzK+UFgvoDX3CCSosQd/6m2/q9ZpNWArNZWluVvddLjcYqBuua9fSCUumHJZSkW1k2F136TmdODEd10oizvLMM/pA06dP7qSQGTNURIuJaTU16WscQYjKYfEUIyyjR4esy7Dl4Ed/BeMmW7cy5eD0TbBXr2aG7ZuHxRJXWOL6dd99F774xfY3ZkuLilrw9/iclR494u03AapGWESkh4gUO5Xjh4AzgIOBC0UkZgnA0pCvxeKfysaMyW2x5BpqnCSzZmlnEJzyOFfgHvQJvKEh/05/3jwduRNHWOrr1bovlSusUsIS7rjXr1exPe007RDHj9d+qth4x65d8Oyz6sqLw4wZ6nr6xz8KP6avGBFHEKJyWDz77KMPY4UIyx6WcNhiySAsXSUtvscctYN6QhZL3OD9E0/oa5SwRAUl165Nl95/7TW9CGbOhF/+UpU+uK0JC4hInYicDbwFTA99N0xEbhCRi0TkdhGZlGN3dzrndjrnNgFvABXIE04zeLCOWoyTJb1kSToIPmaMPoxkq7OVLTkyaT72MRWwYBB//nx9KPNziUQhUtjIsDlz9D738aZc+JFhpcA/FJRTWDZsSHe+nvvv1wfT007Tz/X1WousWGF5/XXtd448Mt76ScRZgsKSy/LJZrFAYQH8WMLiO+OgWOzapU88Kc45c7vepHV1aYslrivMX9xRWfdRwvLd76qvde7cdIE6H5gNBtpaWtq7wrp27ZzCAgwEngTaXToiIsBfgLudczcAPwLuE5Eue+5Ccc61pLYdBPzVOZel8lvp6dJFxSWXxeIDlP6pbOxYdTOFO5cg5bRYfBD/oYfSN/orr+goqGCMMIrx4/Pv9OfM0WTOuK7gCRM0QbEUw2CXLdP7MVUdveT40i5hd5hPaJ02Lb1s6lT9PxTjlnrmGX2Na7GMG6fXdDFxFn9dv/9+7nCCT04NJ2x6Jk/WcxV3EIPPEdtDWMKuMN/IYNwE2pU7+Mi/7lBhiLJY4gbv41osvory73+/Z0JVcP2wxdK3b7rNnUlYnHOrnXPvRnx1LDARaEyt9wbQCpwiIseLyF8j/vqlBOkk4Mdl+glZiZPLsnq1PuwELRbI7g4rp8UCKiyQLqE+f372+IpnwgTtHOLWmFqzRoe+xnGDefbfX++nUlQG9jksxYz2yoeokWFbt2pxy1NOad+nTJminWQho+48c+eqUIweHW99EbVaihGWpqZ0H5fL2li6NHNyKqiwtLbGd81FjgiDPS0WjxcL3+Bg5vKOgLAUGmOJa7EMHKivr7yy58UYZbH4Exa0WKLaVCIqLixZOApY7K2QFG8BH3bOzXHOHRvxtwk4Ffijc26XiJTJgZGZOMLirYBgjAVyC0tDQ+YqEUkzapTWF7vlFu3I1q3LHl/xTJigQhS3SoUv45KPsCRRlywT2ZIjS8Hgwfp/DQrLww9rH+bdYB4v7MUE8J95Jl2uJS4zZug1W6iQr1oFH/ygvs8lLJlyWDz+GozrEowlLEHzyAuLtxiCBdp2RFgs4cz7oEssbvA+Sli8Wbp1a3aLxQfvzzhDP++7b0VcYRndSlXAYCA8sH4DkDGjQEQuAr4NbBaRbsAlwLIM684CZgHsu+++NDY2FtTILVu25Nh2HEuWDKax8W8Z13jyyUHAQaxe/QKNjVtpaakDPshjjy1m+PDou/fVVw+if/9eNDa+UFC7C+HIIwfy0EOT+NKXVgJDce5lGhuzD1HduLEPMI277nqdtTESDm6/fX/69h3I5s3P4E9rrnO8fn1X4CgeeGARPXvmUZwtBosWHcGhh66nsbGIaHWeDBs2lXnzHI2N2lvecMOB9Ou3F21tc2lsTPu9mpvrqKubwT33LGGvvd7N+zhr13bj3XeP5GMfW0Rj4/IY17LSs6f+T2+6aSHHHZd/zsXy5UczceIqBg/em0cf3cwRR2TOCP3nPw9n8uQNGc//rl1Ct25Hc//9Kxk5MvfTy1//OgLYjxUr/sbGjWkX1/5r1zIk9X5b9+54CVj+zjsMaWmhuUsXegFbV67E58rPnzePwcuXs9euXWzYsIG+Gzaw/O9/ZwLwblMTo4E253Y/vb+1bBmpZ6Dd5/no1tY9OuHtmzfTM7TeiEWL2A/YtnYta5YvJ5iTvHDBAt5LieGkpiZ6NDfz4rHH0vVf/oUDr7qKvVLC8upbb7GuwH4ub5xzVfEHOODYwOfrgKdC6/wO+EvSxz700ENdoTzxxBNZv7/6aufAuY0bM69zzTW6zrp16WVDhzp3zjmZtznsMOeOOy6/thZLS4tzQ4Y4J6J/mzfn3mbjRv1tP/pR7nXb2pwbPty508EXO8AAACAASURBVE9vvzzXOW5rc66hwbmLLsp9jHxobXWurs6573432f3m4gtfcG7QIH3f3Oxcv37OnXde9LoHHeTciScWdpw//Un/N88+q59znWfPzp3aplmz8j/mtm16zKuucu7kk5074IDM67a2Oldf79x3vpN9n9OmOXfMMfGOf9FFeq3swfnna8PAuUMOSb+/6CLnunXTGw6cGzEi/d2cOXqTjhihryNHOnfttfrdlVfqa7du6fVvvTX93tOrV3qZ/xs2bM/1fvzj9HeXXdZ+/dmz0+t99KPOTZ+e/nziien1Hn883kmKCfCiy9CnVrMrrAkIDz9uAIrIHio/ceZlWbJEY2zB+dFzDTmOmpK41HTtmi5fMn68TnaVi3791BqP46b6xz/UbZiPGwzSo8+SHnLc1KSeimxZ96Vg4kSNNa1Zo8UhN23KPL/LlCmFjwx75hn1jsSJlQWpr9eSMoXEWXxMfMgQjY+89Vbm+FtTk3qacrki/ciwOIMYIkeEQXtXWNBl5Odg8SNJsrnC8i1CCdGusKjZJoOusLDfMip47wm+z5VtnCDVLCxPAGNFJNjGA0gF82sFH1zPFmfxI8KC10u2JMldu/SmK1fgPsgXvqDtjBNf8UyYEE9Y4pZxKeYY+VDuHBaPD+AvXKijwfr0gWOPjV536lS9tgqZumPuXK3uke/cMaBxloUL86/R5oVl8GAVhLY2HfIcRbYcliCTJ+vv94mX2chbWPyskf1Tz7iZhKV79+g8FhEtVxBcFiRKWKIqDmeLsUQF7z3+/ciR8UdoJEBVCEtIPDzzgCXAzNQ6BwC9gRiFIGIf9yQRuWljnBr1BRInSTKYw+IZM0a3iRpGuWaNXs/ltlhAr83Zs3UWx7jEzWWZM0dHeGUqapmN/ffX8xWcJrhYKi0sr72mc+d8/OOZ4675Bq8927bp6Lu4+SthfD7L3zKHDiMJCwtkDuDnymHxxM3Ab2vTgSeRwhIcbhw82V5IvMXiXPr7KIslalTYvHnQ2BhdFyxqVFg2i6W1dU+LJWjyZbJYiqmOWgAVFxYR2QcNuAOcJSIHgjoXgZOBs0XkEuAy4ETnXGJdh3PuPufcrP79i034z4zv/ONYLEHGjEnnt4Qp91DjMOeck5/FMn681tzKNuNhc7Pee/m6wTzF1CXLRKWEZfhwtVJ+9St9iAiPBgtSaGmXF15Qr03c/JUw06Zp35mvO8xbFUOG6DXep09mQYh7/r2w5BLXVav0OsvLYvGmYDCpyvuAoyyWnTu14w9aDYMG6VzfUcKSa1TY+tQkd0E/X/hpMygsYYvF78snSJWJiguL0zyWHzrnxDl3rnNuYeC7t51zZzvnrk+9lm8IVEL06KHXVaYYy9atmigWFpaxY/U1yh1WaWHJlzid/ty5+hRdiBsMSlOMcvly7UNK+NwRiYhaLQsWaH91wgmZ1917b+1487VY5s7V1yOOKKyN3bvD4YfnLyyrVqknZ9AgfT3kkOwWS79+uc9///5qSeeyWDIONYbMwuLn9R40KN1h+zlUtm9vb7Hs3KnK1aVLdOJN1LJcgaG99lJfYXC9sFkeFpagxeKfaL07rkxUXFg6A9lyWbxFEuUKg+gAfjmz7pMgTjHKfMu4hBk3Lvcx8sXPHFmu5Mgg3h123HG5J2ebMiV/i+WZZ7QcT6bZR+MwY4Yed3PmCT73YNUq7aP9w/vkyVolO6p/zSeHKE5pl6KEpVu3dIwiymLxnfnWrSpA/qIJXjyZsjxBa35l4tZb25+gQFkZYE9XWNBi8R2MCUv5KEeMBdSyyCQsPsksbLEMHarXR5SwrFih12i7st9VzH776f2VS1iOPDK/GS6D9O6tIpDkyLByJ0cG8Z6LbG4wz9Sp+rvD/U0m2trUYinUDeaZMUP71Xnz4m+zapXGVzyTJ2sye1SyZZSLOBNHHaWjCrNdY/5eioxhZ4qxeFdYfX3ajeAtlqCw+O2bm9tnvgeFJcoV5sk2xHLduvjCErZYvD/RhKV8lCPGAoVZLPX1egNEucJWrlRRyVSFotro2VM7iEyusCVL8i/jEkXSxSgrKSwnnaTnI9Mw4yBTpqhYvPZavH2/+aa67osVliOO0H7z2Wfjb9PUpPEVT7bAez7n/3Of0/vhV7/KvM7ixSpqUSW7MlosnqCwdOum4hEUFr99c7N+jusK82QTlvXr47vCwsH7O+/UIYXlKnaXolMLS7kYPlwt6qjx+kuW6HUYvNk8mXJZ4kxJXG1EFaPctQuuu047l27d4j2dZ8PnshRTlNHT0qIDDippsTzySPvcpkz4PJS4cRZfeLLQEWGevn31uo0z26knbLEcfLCKU1hYtm1TYyGuxTJ4MHziE3DbbdGDqiDLUGNoLyxRw4KDwrJ5s4pPlMXiZ3CMEpGkLJZcMZagK+zUU+HRRzPvu0SYsJQBP+Q4qpKrL7IXZX1kEpZKJEcWi88z8ffHc89pDsWXvgTTp6ufvdiBKxMmqFsl39yKKFau1LZWSljyYdQoFaB8hGXvvdODKoo9dtyaYW1tKtZBYendW+NjYWEpZLqCWbNUjO69N/r7rMIS7PR9xxx2Y/kL9M031ex54w0Vlrq69A3c2qqfo1xh2SyWbD7gfCyWsCusQpiwlIFsuSxLlmR+Khs7Vi2d8DDdWrRYfKf/5ptw4YXqRlm1Cv7wh3T+SrEkOTLMu6bLnXVfCD5hNW4Af+5ctVaSGJQwenT86srr1mm/G7bOowLv4cKscTjuOBW64LxBntZW/Z/Gsli8sARHNtTXw0EH6fsdO7Tc9yOPqB8wbLEEhSVINmHxcZsowhZLeH6IBQvg2mvTx881l0UZ6NTCUs7gPUQLy9KlmRMCo0aGbd+u11mtCYsfGTZ1qlZI/trXNNh65pnJjbryT+BJBPArlcNSKFOmqNUXlW8XZM0aFd5i4yueUaP0XOU6LrRPjgwyebJWv96yJb2skPNfV6eVIR57bM9q2suWqcWUl7AE/ZD19e1V7pRT9HXbtvjCks0VFhn4SRG2WMLVj5cvh3//d30C3batLHPa56JTC0u5gveZhGXXLl2W6anM3wTBAL5PMKs1V9ghh6iFfuihGqi/5prCR4BlYvRo7ROStFhqSVi2b8/9233+SlLCMnq0PkDHKaeSTVig/eADb7HkazGee6723zff3H551qHGEC0sQZeSj5v07681/4MXbzB4X6jF4svyR5nu27fDn/6U/pxpRjvfwcQJzJWYTi0s5aJvX70ew8KycqWKSyaLxccKgxZLrSVHeoYP187nqadUZEpBfb3665MSlv79kxe/UhE3gD93rvabhx6azHH9tRvHHZZLWILusGXLdORjvnNTDRumJXBmz26foJ5TWKJiLEEh8N+vW6clIsLC4r9vbW0fvI873Lh3b7VKrrgi+vvgTGaZpsv0alzujN4ITFjKxPDhe2bf5/IjDxigVm1QWGotOTLIXntlf2hLgqSqHC9fXjvWCmhcuVu33HGWZ55RUcnmeckHLyxxAvjBci5BRozQh+ygsOSTwxJm1iydlfW+QFXBxYu1X8/4P42yWKKEpa5OxSKTxRIO3gfJdvH77ePER6KKVELazDZh6TxE5bL4mzGTxSKiT1hBV1itWizlYv/9YdGieD7/bPis+1qha1eYNCm7xdLcDC++mJwbDPITllWrtHJ8eGStyJ6lXYrJIfroR/V/FwziL16s+8uY+xVXWDzBH5FE8N5bNuERXYcdtmfdnUzC4p9UzRXWeYjKvo8z8mXs2D1dYT16VMW1U5VMmKAPjcXMAw+VTY4sFD8yLFMez8svq7gUm78SpFcvLdES1xU2eHD0YA1f2qWtLV18tVCLpb5eB209+mj63sk61Nhv5PEik01Y6uvT86vEtViyucI8YYulvl5PcJBcwmIWS2Up16gw0CeoVavaD+hYskRHNGYbaehzWXxnsXKlilQl6lfVAn5kWDFxlh07dPRUrQnL1KmaxxGVLwXpxMgkLRbQAH5ciyUqERhUWLZuVet8/Xod3FTM+T/vPL1HbrlFP+cUlnwtFki7w6IslqgbNEpszjuv/ecoYQlbM2axVDflGhUGKizOtR89E+epbOxYvY584LMWc1jKSRK5LIUk51UDmUrob9gAP/gBXHWVnp+ka8yNGhXPYmlq2jNw7wkG8AvJYQkzcqRWhZ49W0fhvvdeAcISFJO4wpItQTJqH7fc0t7EzCYsIukqylFYjKXzETVFcdQEX2HCuSy1mHVfTgYN0vuqmAB+rQpLeF6StWvhO9/Ra+yKK3SU7J13Jn/cUaNUDHKV0gmXcwly0EHaFy9YkNxQ7wsuUDG7/nr9XFKLpdjM+/Cxg9uEhSUT5grrfISz753LnnXvCQqLc2ax5EKk+GKUtZR1H6RvXx1u3dgIl16qLqr//E8tZvnKK1rqxCePJ8no0WpVr16deZ3mZnVxZRKWnj31/5aUxQI67HjIELj6av0cO8YSHAEW9b0nmyss3+C9J47FkonWVvWrV0F1WhOWMhFOktywQTONc1ksvsT3O+/oNjt2mMWSC1+XrFBqLTkyyJQpmnn+059qcvjrr2tunbdmSkGcXBbvys0UY4F0aZdly7R/LdZl16WLBvE3bNDPWYXFd+j9+8d3Y0VZLIVm3ofbEdwmm7D4AQSeKrBWwISlbAwYoE9lXljiPpX17JmuIGtDjeMxYYJ2TuFafXFZtkxzbsL3bC3w1a9qdY9//AN++9v0hGGlxD/8ZAvgZ0qODDJ5su7jtdfUWkwi5+n887U/7tEj+7GZORO+/31VtrgWiy+dkq26caldYWHrpAoC99DJhaWco8JE2uey5MphCeJHhnlhMYslOz6Av2hRYdvX4lBjz1FHabkcX5utHORjseQSFoDHH0/u/I8erS6xSZNyjKTs2RMuv1x/jBeAXMLirYNswfsgpbBYwmadWSyVp5yjwqB99n0+fuSxY9UV5oeRmsWSnWKLUdZa1n2l6ddPH5TjWCy5XGGg7t5i4ytBfvc7ePjhPDbwopBrVJi3DoJl83fuTD7GErR+wsJy8MFaJ8nP6W3C0vkIWyzdu8M+++TebswY3c4/EZrFkh3/tF5onKXWsu6rgVy5LE1N2i+Gc/2CDBkCAwfq+ySFvW/f9hXwcxLXYvHC4meNDG5fqLCEhSOXK6y+XueI9qagCUvnY9gwtVja2tI5LHESHceM0W2efVZ9/1EzpxppevdWYShEWLZt0zqDZrHkR65cllWrVDSylcISSVstSVoseRMlAFEjrbywbN7c/vt88ljCRAlHnBiLn3rYhKXzMXy4WsqrV8fLYfH4Ksfz5pkbLC6FFqOs5RFhlcRbLJlyWbLlsATxwlLR8+9Foa0tvSybxbJ5c3IWS64YS9T3kC67XwVzsYAJS1kJ5rLkUwvJD5PcvNncYHEZN6598c64mLAUxqhROnx+3bro75uassdXPIcdpq/lHHywB5UUlnxdYd5i8XWhSl0+PCbV0YpOgheWd97RGy2usAwbln5QMYslHgMH7jnxXhxqNeu+0uSqchzXYjn9dK3A7AdgVATfOQdLZOcSlkyusCClGBXmj+v3HRTDCmLCUka8sDz7rL7GdYXV16fXNWGJR0OD9gvB6W7j4C0WO8/5kS2Xxbn4wlJXl9wkZAXjVe3zn08vy9diCYpBcHkuwrGcXMISzrkpdr6IhOjUwlLOPBbQETFdu6arzOYToPTuMHOFxcPf8z7rOi7Llun/yQZI5Ee2XJYNGzRvMI4rrCoYMkTV8Pzz08sKtVjyDd6HR/OEKyVnslj8cDs/rK7CdGphKXceS12dCsPLL+vnuBYLpIXFnqTj4WOZhQiLucHyZ6+9dO6rKIvFV/SOY7FULdmEpUePPXNeCo2xhBGJZ7F89rNw443wjW/kf4wS0KmFpRL4kWH+fVz8yDCzWOLh7/n16/PbrqnJxLsQRDIPOY6TdV/1RAlL9+5w3XVa9TNOjKVQYcmWIBmMsVx4YfYilWXEhKXMeDEZMkSvy7h85COaXHvAASVpVoejUFfYunXplAAjP0aNirZYOqywAFxyicZkkhoVFiaXxVIFlYyjMGEpM15Y8k0AmzIFnngi+2yTRppChWX9+rQbzciPTNn3ccq5VD254iOZhCUYH0m9d/lO/xrHFVZlVKfcdWC8myWf+IqRP4UIS2urjiIzYSmMUaNUmDdtap+n19SkYYgqyd0rjFwdeBxXGMDVV/PiXnsxPe5xw8LUES0WERkoIt8QkXkislpEtovIIhG5Q0Q+VqpGdiQKtViM/ChEWPy6JiyFkWnIsR9qnO+DelVRqMUS5tJL2brffvGPG3aFhfNc3nsv/r7KSGy5E5GzgbOA54H/BjYBO4B+wGDg06l1LnbOvV+CtnYIvLCYxVJaunTRUUr5BO/9uiYshRFMkjz44PTyVatq3A0GyVkshRAlLCNHwnHHwTnnJHecBIklLCJyEbDEOXd8ltVuEJEhwDdF5HvOue2JtLCDccghOj/E8dnOpJEIDQ35WSwmLMWRKZelqanCmfRJkESMpRDCrjAvYN27w803F7fvEpJTVkVkb+Bp59yDudZ1zjUB3wEOzrVuZ6V3b7j//g5wo9UAJizlZd99NZaSyRVW0yTlCsuXTK6wKindkomcv945975z7nWIF2Nxzu10zj1f6oYnQbkz743yUqiwVMnsrjWHiHpoghZLSwu8/34nEJZghnymBMlCiRKWKindkonYvz4VP/kd0IDGWM4FTgQuBZ5AYyx/TFk4NUG5M++N8jJggFks5SY85Hj1an3t8DGW4DrhMizFEE6QrBFhsRiL0WFpaIBXX42/vglL8YwaBa+8kv7cIcq5QDxh6dJlz6mJkxCYjmixWIzFqFUKcYX16GEFKIth1Ci1UrZt088dIuse4rm2oiyWpIXFB+87UozFIyK3iEhkcetairEYHZuGBk3Wi3sPWtZ98fhclqVL9bVDZN3HJVzCPgk6avA+A4OBpiQbYhhJM2CAVj/ftCne+iYsxRMecuyFZZ99KtKc8uKtiXIIS627wjJwKXBCKqZiGFVJvtn3GzaYsBRLOPu+qQn23rtqiu6WlqDFku/UpZnoZMLyACouT4rIz0TkoATbZBiJkG/pfLNYimfIEH1w98LSIXJY4hJlsVjwPi++ARwLXAnsCzwlIgsTa5VhJEC+FosJS/HU1+tEaUFXWKeIr0B7i8VffD/4QXH7rNEYS0GlMZ1zd6be/gb4jYjUAYcn1irDSAATlsoQzGVpaoKjj65oc8pH0GLp3j0Zd1hYWPwxat1iEZE+ItI12zrOuTbn3NzANnZ7GhUnn+mJd+3SIL8JS/H4mSSd6wCusL3zyPf2FkvSc6R0UFfYDuAyEekbZ4cicg4wrphGGUYS5BNjsZL5yTF6tFoqa9fCjh017gp74434WbalGhVWg5n3cfJYdgK/AH4tIheLyNjwOqkaYieLyN3ABufcCyVoq2HkRb9+ei/GsVgs6z45Ro1Sa+WFVC9Q0xbLPvu0nwMgG0kJyxNPpN/XaIwl1hlwzq0FzgAGAI0i0pwqQtkkItuAt4GTga845+4pXXMNIz51dSouJizlxeeyPPusvta0sOSD7/SLFZaZM+Gaa9Kfs030VaXEDt4753aKyK+BHwKHACOB7sAy4DXn3LbSNNEwCiduIUoTluTwuSwmLEXgA/+ZgvdVTj4zSArwOjAC+BFwsnOupVQNKwcichJw0rhxFhLqqMStF2bCkhzDh2vf+nyqsFNNx1jyIUlh8eSamrhKyecM9ACeBiYAk4EDgl+KyJUi8gcRuUxEBibYxpJhZfM7Pg0N8YL3JizJ0bUrDB0KGzdqxn2nmd/GlxdI0mKBmhSWfOyqk4GngA8A9wJrQ99/ERgCtAH3AR9PooGGUQwNDbBoUe71TFiSZfRoWL5c3WBJTU1S9ZTDFVYjwpLPGbgf2AT8A3gNFRpgt5usX2oE2YHAtCQbaRiFko8rrHt36Nmz9G3qDPgAfqdxg0Gyw41rXFjyCd5vAX6Z+vi0iEwXkR8AdwP/CqxKrfd3qx1mVAv5BO/NWkkOH8DvNIF7SAtKKRMkO1rwPoxz7gUReQv4LFpG/7TAd2E3mWFUhIYG2LJFJ/bLdk+asCSLt1g6lbAkOR9L0GIJ7q8jWSwi0g/Y7pxrTX0eDuzvnHsMuL6E7TOMogjWCxuYZUiJCUuymLAUSZQrDGpGWOLUCvse8CbwkIh8W0TEObccExSjBohbiNKEJVn2209fhw+vbDvKSqlnkISaEZY4FsvZwATn3ObUMOKzReQOoLW0TTOM4olbiHL9ejjIIoOJsd9+cM89cOyxlW5JGUkyeN+nj7727dt+6HEHEpaXga2wO3Zym4icCvQpZcMMIwnMYqkcJ5+ce50ORZIWy4UXwvbt8OUvw09+kl7egYL3FwOni8idqeHEOOfuFpF9S9s0wyieOMKya5cm85mwGEWRpLB07QqXXqrvg66wpEeclYicwuKcew/4Y8TyG0vSIsNIkDil8zdu1FcTFqMoShFjgZrMMC3oDIjILSJyaNKNMYykiWOxWNa9kQgmLLsp9AwMBpqSbIhhlII+ffR+N2ExSo4XFBOWgoXlUuAEEelMBRuMGkQkd1kXExYjEbzFksRc90GihGXo0GSPkTCFDjF4AGgGviUiDwK/cs79PblmGUZy5KpwbMJiJIIfsbVzZ7L79RaQF6w33oBBg5I9RsIUKizfAJ4DPgR8DHhKRFY75w5MrGWGkRBmsRhlwVssSc9HH7ZYDqz+bjansIhIH6DZl3MBcM7dmXr7G+A3IlIHHB7YZoBzLsYsGIZRenIVojRhMRKhXMJSA8SJsewALhORvplWcM61OefmAojIOUBFpmQUkT4icq2I/FVEvlWJNhjVRxyLpVs3K5lvFIkJy25yCksqKfIXwK9F5GIRGRteR0QGisjJInI3sME590K+DRGRHiJS7FSOY1E33fHAcUXuy+ggxBGWAQNq8v41qgkTlt3EGhWWKuVyBjAAaBSRZhFZLSJNIrINeBud+Osrzrl78mmAiNSJyNnAW8D00HfDROQGEblIRG4XkUk52vmqc24X6pb7VT7tMDoucYL35gYzisaEZTf5TPS1E7hKRH4IHAKMBLoDy4DXnHPbCmzDQOBJYERwYWpWyr8Alznn5ojIk8ADIjLel5aJQkRGAhcCh4nIvc65HQW2y+ggNDTAjh3616PHnt+bsBiJYMKym7xHhTnnHLAg9bcbETnSx1ny3N/q1Pbhr44FJgKNqfXeEJFW4BQR2QR8M2J3pznnlqIVmH8PHAzk7ZYzOhZeNDZuzCwsnWreEKM0mLDspqhSmSLSFfg08GVgarH7C3EUsNg51xJY9hbwYefcxcCcHNs3Ae8k2B6jRgmWddk3onTq+vU1MYLTqHZKLSxJJ16WkIKEIFXZ+CLU5TQIWAIkLauDgY2hZRuAjFMHicjXUYH7A/Cgc+79LOvOAmYB7LvvvjQ2NhbUyC1bthS8rRGPYs/x0qV7AYfw2GMv0dS0eY/v16w5im3b3qOxcVHhjewA2LVcHCOWLGE/YOnixbyT5Tzme56HLlrEhNR2L9bI/ycvYRGR6cBXgE+hnf4twI1okmTSM0ruZM/JxLIONnDOXRN35865m4CbAKZNm+ZmzpyZb/sAaGxspNBtjXgUe467ddPXsWMPJbybtjbYuhUmTRrOzJmdabrDPbFruUhefBGAkcOGMTLLecz7PC9cCECf3r1r5v8Ta1SYiHxGROYBz6JDes8HRjjn/m8qplEKG60JCA8/bgBWluBYRgcmW4XjTZvUw2DBe6NoLMaym3yKUO5C52U50Tn3v6HYRyl4Ahibyur3HEAqmG8Ycck2PbFl3RuJYcKym7h5LL93zh0N/Bj4nohcl+R8LCHx8MxDYzczU+scAPQG7kvwuCeJyE0bN4ZDOUZHIpvFYsJiJMYnPqGv552X7H5rUFjyirE45xYAXxaRBuBcEfkK8Aiaz1IQIrIP8IXUx7NEZIVzbqFzzonIycDlInIgcBhqLW0t9FhhnHP3AfdNmzbtgqT2aVQfPXponCUqSdIv8+JjGAUzenRpRm51dGHxOOc2ANcCiMhH0GB+QWc0lcfyw9Rf+Lu3gbNTH5MeHGB0ErLNyWIWi1H1dBZhCeKcewR4RERGJdAewygJmSocm7AYVU8NCktic2g655Ykta9yYTGWzoNZLEbN0pmFpRZxzt3nnJvVv3+xRZWNaiebsHTpAr17l79NhhGLGsy879TCYnQeMlU4tpL5RtVTV3vddO212DAKIJvFYm4wo6qpwaceExajU+CD92FvggmLUfWYsNQWFrzvPDQ0QGsrbN/efrkJi1H1mLDUFha87zxkyr43YTGqHhMWw6hOvLCEA/gmLEbVY8JiGNVJVCHKtjb9bMJiVDUmLIZRnUS5wjZvVnExYTGqGhOW2sKC952HKGGxrHujJjBhqS0seN95MGExahafIGmZ94ZRXUQF701YjJrALBbDqE66dYNevcxiMWoQExbDqF7CZV1MWIyawITFMKoXExajJjFhMYzqJVzheP16qK+HPn0q1ybDyIkJS21hw407F+FZJK1kvlET1OAF2qmFxYYbdy6iXGHmBjOqHhMWw6heTFiMmsSExTCqFy8sPs/MhMWoCSxB0jCql4YGrQ22ebN+NmExagKzWAyjeglXODZhMWoCExbDqF6C9cKcs5L5Ro1gwmIY1UtQWDZvhl27TFiMGsCEpbawPJbORVBYLOveqBlMWGoLy2PpXHgRWb/ehMWoIUxYDKN6MYvFqElMWAyjeunXT19NWIyawoTFMKqXLl2gb18VFj/k2ITFqHpMWAyjuvEVjs1iMWoGy7w3jOrGVzj2JfP79q10iwwjB2axGEZ14+uFrV+v72vwnjU6GzV4kZqwGJ2KoLCYG8yoCUxYDKO6MWExag4TltrCMu87+nNh5gAACopJREFUH8HgvQmLUROYsNQWlnnf+RgwADZtgrVrTViMGsGExTCqG599v2yZCYtRI5iwGEZ144WlpcWExagRTFgMo7rxwgImLEaNYAmShlHdBMXEhMWoCbzFYsJiGNWJWSxGzWGuMMOoboLCEnxvGFWLCYthVDdmsRg1hwmLYVQ3ffum71MTFqMmMGExjOqmri5ttZiwGDWBCYthVD++qrGfUdIwqhoTFsOofhoa9K/Orn6jFjBhMYzqp6HB3GBGDVGDwtKl0g0wjHJz5JEwcmSlW2EYManBzPtOLSwichJw0rhx4yrdFKOM/Md/VLoFhpEHNWixdGpXmJXNNwyj6jFhMQzDMBLFhMUwDMNIFBMWwzAMo7NjwmIYhmEkigmLYRiGkSgmLIZhGEaimLAYhmHUAjWUIGnCYhiGUc3YqDDDMAyjs2PCYhiGYSSKCYthGIaRKCYshmEYRqKYsBiGYRiJYsJiGIZhJIoJi2EYhpEoJiyGYRi1gCVIVhYRmSoiN1a6HYZhGJ2RqhEWEekhIkVP5SgifYEPAz2Kb5VhGEaFscz7/BGROhE5G3gLmB76bpiI3CAiF4nI7SIyKcYuPwncVYq2GoZhGLnpUukGAAOBJ4ERwYUiIsBfgMucc3NE5EngAREZ75zbGbUjETkReAjoWeI2G4ZhlIe+ffV12rTKtiMPKi4szrnVALKnuXcsMBFoTK33hoi0AqeIyCbgmxG72wx8AegFHCAiX3XO/axETTcMwyg9gwfD88/DpDgOm+qg4sKShaOAxc65lsCyt4APO+cuBuZk2lBERgPfM1ExDKNDMH167nWqiGoWlsHAxtCyDcDwJHYuIrOAWQD77rsvjY2NBe1ny5YtBW9rxMPOcXmw81weOsN5rmZh2Qm0hpbFGmzgnHsXOCfHOjcBNwFMmzbNzZw5M+8GAjQ2NlLotkY87ByXBzvP5aEznOeKjwrLQhMQHn7cAKysQFsMwzCMmFSzsDwBjBWRYBsPIBXMNwzDMKqTqhCWkHh45gFLgJmpdQ4AegP3JXjck0Tkpo0bw6EcwzAMo1AqLiwisg/w7dTHs0TkQADnnANOBs4WkUuAy4ATnXNbkzq2c+4+59ys/v2LTvg3DMMwUlQ8eJ/KY/lh6i/83dvA2amP15ezXYZhGEZhVNxiMQzDMDoW4mqoFHOpEJE1aI5MONjSP8aygcDa0rUuZ3tKtX2udQv9Ps45tXMcb51iznHUMjvP+X1XyLUM5TvPpT7Ho5xzgyK/cc7Zn4rrTYUsA16sZBtLtX2udQv9PuY5tXMcY51izrGd5/Kc5wzrlOU8l/Mch//MFZYmarRZ3GXlothj57N9rnUL/T7OObVzHG+dYs5x3OOXis5ynjvLOW6HucKKRERedM7VTtnRGsTOcXmw81weOsN5NouleG6qdAM6AXaOy4Od5/LQ4c+zWSyGYRhGopjFUgZE5OBKt8EwisGuYSMfTFhKjIgcBjxb6XZ0NESkq4h8T0ROFZFLK92ejoxdw6VFRAaLyF0iskREvl/p9iSBCUuJcc49D6ypdDtqARHpISJx6+t8AXjbOXc30EdEji5h0zo1dg3nT57X8oeAM4CDgQtFpKF0LSsPJiwZyPPCMIpAROpE5Gx0htDpoe+GicgNInKRiNwuIn5+1sOBV1PvXwE+Wr4W1zZ2bZeOAq/lO51zO51zm4A3gO1lbnbimLCEKPDCMIpjIPAkMCK4UEQE+Atwt3PuBuBHwH0i0gWdYXRLatXNQHQGsLEbu7bLQt7XsktNvy4ig4C/Oueay9zmxKl4EcoqJNeFcZlzbo6IPAk8ICLjgQ8D34zY12mppxAjC04LkaKnuB3HAhNJzcHjnHtDRFqBU4D3gT6p9fpQvlIktUze17ZzbmcF2lmzFHgt/zn1PzgJ+HHZGltCTFhCFHJhOOf+DMwpYzM7C0cBi/0TXYq3UCF/AvVJLwAOAR4rf/Nqi0I7vTI2sSOT7Vr+M3Aq8Efn3C4RGeGcW1aJRiaFucLik+3CyIiIfAAYJCLHl7JxHZTB7FkEbwMwHLgVOFBEzkCn73m83I3rQGS9tu0aToSM17KIXARcCzwnIm+hM+XWNGaxxCdbJ5cR59zL6MyXRv7sBFpDy+oAUi6a75S9RR2TrNe2XcOJkO1avgG4oewtKiFmscQn44VhlIwmtHR3kAZgZQXa0pGxa7v0dKpr2S6e+HSqC6NKeAIYKyLB6/QAUrEAIzHs2i49nepaNmGJT6e6MMpN6Lx65gFLgJmpdQ5AXTKVLEXeEbFrO0HsWrYYSyQxLozHO/qFUU5EZB80kx7gLBFZ4Zxb6JxzInIycLmIHAgcBpzonNtascbWOHZtlxa7lhWrbhwicGFcBdwGXO2cW5j6bj/gcuB59MK4zjn3QoWaahh5Yde2US5MWAzDMIxEsRiLYRiGkSgmLIZhGEaimLAYhmEYiWLCYhiGYSSKCYthGIaRKCYshmEYRqKYsBiGYRiJYsJiGIZhJIoJi2EYhpEoJiyGUWZEpKuI9C3j8QaW61iGASYshlFWRGQQOkFZi4h8XEReFxEnIhMj1u0pIutEZKuInFfEYYeIyMVFbG8YeWHCYhhlQkS6AjcB1zrnmp1zDwD3ohNt/Z+ITc5MfTffOTe70OM6514DlorIOYXuwzDywYTFMMrHt4DHnHPBaYBbgbuBz4lIv9D6JwKPoOJSFM65+9Ey7vsUuy/DyIUJi2GUgdQ8KBcCD0d8/QugF3BOYP0pwKvArtB+zhORF0XkHBFZICLvicgFwe1E5CoR+baI/FVExgY2fx44O7lfZRjRmLAYRkKIyKdE5F0ROUhEHhSRWwNfHwIMds4titj0XeAB4BIRkdSy84GbI9a9DzgUnfJiMvAj4AYRGScig4EbgCuccz8C3kfFzPN34FNF/ETDiIXNIGkYyfEc0AwcCcwChgS+GwlsjNooxXWo2+t4EXkG6OOcW5nWGcU5tya17InUov8BvgscB/QFnnPOedfZ2UBwwqVNwH75/yzDyA8TFsNIjg8B9cCDzrkVwPLAd90JubVCPAq8iQbxxwC/jXNA59xOEVkM9ERFY1vgux2h1Xeg4mMYJcVcYYaRHB8jLSphlgH9M23odCrXX6T2cQrweB7H7QMsBFaj89bvRkTGBT72A6LaZhiJYsJiGAkgIvXA8UTHRQBeA1ojRmX1RoUBdB76bcAjLj1neDega8T+eqeOOwEdNTYH+DMwWUR+LiITROR0YEJgm8HAy/n8LsMoBBMWw0iGw4E1zrlXo750zm0F7gSO8stE5GTgM8BlIjLOObcJ+CUqMIjIJ1ELZIqInBva5RdF5HvAFcAnnXO7nHML0KD/qcDfgJHOuQcD2xzh920YpUTSD0aGYZQSERkF/Ldz7uQi9+OAMc65d/PYZijwc+AMZze9UWLMYjGMMuGcWwL8usjyLB7JvUpqRc34/xpwsYmKUQ5MWAyjjDjn7gL+LiJ7F7K9iHwx9facVN5KHD4AfM85t6aQYxpGvpgrzDAMw0gUs1gMwzCMRDFhMQzDMBLFhMUwDMNIFBMWwzAMI1FMWAzDMIxEMWExDMMwEsWExTAMw0gUExbDMAwjUf4/SsY0M7oHkeIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "max diff for r = 0.1-100:  0.02554546783902975\n",
      "max diff for r = 50-250:  0.027782293661027104\n",
      "average Delta(r^2 xi) for r=140-150 Mpc: 0.006977835181869886\n"
     ]
    }
   ],
   "source": [
    "# plot absolute difference in r^2 xi(r)\n",
    "plt.plot(r1, abs_diff1, 'b-')\n",
    "plt.plot(r2, abs_diff2, 'r-')\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.xlabel(r'$r$ (Mpc)')\n",
    "plt.ylabel(r'$\\Delta (r^2 \\xi(r)) $')\n",
    "plt.grid(which='minor')\n",
    "plt.title('Absolute difference') \n",
    "plt.grid(which='both')\n",
    "plt.savefig('benchmark_abs.pdf',bbox_inches = 'tight')\n",
    "plt.show()\n",
    "#print abs_diff\n",
    "print(\"max diff for r = 0.1-100: \", np.amax(abs_diff1))\n",
    "print(\"max diff for r = 50-250: \", np.amax(abs_diff2))\n",
    "\n",
    "\n",
    "# find and print the average of Delta(r^2 xi) in the BAO peak region\n",
    "max_value = 0\n",
    "max_value_index = 0\n",
    "avg_value = 0\n",
    "for i in range(63,68):\n",
    "    #print i, r2[i], abs_diff2[i]\n",
    "    avg_value = avg_value + abs_diff2[i]\n",
    "avg_value = avg_value / 5.\n",
    "print(\"average Delta(r^2 xi) for r=140-150 Mpc:\", avg_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for i in range(nr1):\n",
    "    col = []\n",
    "    s = ''\n",
    "    col.append(\"{:.18e}\".format(r1[i]).ljust(25))\n",
    "    col.append(\"{:.18e}\".format(frac_diff1[i]).ljust(25))\n",
    "    col.append(\"{:.18e}\".format(abs_diff1[i]).ljust(25))\n",
    "    s = col[0] + col[1] + col[2]  \n",
    "    #print(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(nr2):\n",
    "    col = []\n",
    "    s = ''\n",
    "    col.append(\"{:.18e}\".format(r2[i]).ljust(25))\n",
    "    col.append(\"{:.18e}\".format(frac_diff2[i]).ljust(25))\n",
    "    col.append(\"{:.18e}\".format(abs_diff2[i]).ljust(25))\n",
    "    s = col[0] + col[1] + col[2]  \n",
    "    #print(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
